

















































































































 


































































































































































 






























































































































































































 



















































 










 








































































































































































































































 










                                
                                

                                
                                

                                
                                

                                

















































































c

      SUBROUTINE ADMmass_integrand3D(origin,Dx,Dy,Dz,x,y,z,gxx,gxy,gxz,
     &     gyy,gyz,gzz,ADMmass_int,Psi,Psi_power,conformal_state)

c
c
c
c
c
c
      
      IMPLICIT NONE

c
      INTEGER,INTENT(IN) ::
     &     Psi_power,conformal_state
      REAL*8,INTENT(IN) :: 
     &     Dx,Dy,Dz,origin(3)
      REAL*8,DIMENSION(:),INTENT(IN) :: 
     &     x,y,z
      REAL*8,DIMENSION(:,:,:),INTENT(IN) :: 
     &     gxx,gxy,gxz,gyy,gyz,gzz,Psi
      
c
      REAL*8,DIMENSION(:,:,:),INTENT(OUT) :: 
     &     ADMmass_int
      
c
      INTEGER :: 
     &     i,j,k,ip
      REAL*8,PARAMETER ::
     &     half = 0.5D0
      REAL*8 :: 
     &     rad,ux,uy,uz,det,dxx,dxy,dxz,dyy,dyz,dzz,uxx,uxy,uxz,uyy,
     &     uyz,uzz,term1,term2,term3,dxx_y,dxx_z,dxy_x,dxy_y,dxy_z,
     &     dyy_x,dxz_x,dxz_y,dxz_z,dyz_x,dzz_x,dyy_z,dyz_y,dyz_z,dzz_y,
     &     Pi,idet,p,pip,pim,pjp,pjm,pkp,pkm

      REAL*8 Cauchy_dr
      integer, parameter :: cctki_use_Cauchy_dr = kind(Cauchy_dr)
      REAL*8 Cauchy_r1
      integer, parameter :: cctki_use_Cauchy_r1 = kind(Cauchy_r1)
      REAL*8 detector1
      integer, parameter :: cctki_use_detector1 = kind(detector1)
      REAL*8 detector2
      integer, parameter :: cctki_use_detector2 = kind(detector2)
      REAL*8 detector3
      integer, parameter :: cctki_use_detector3 = kind(detector3)
      REAL*8 detector4
      integer, parameter :: cctki_use_detector4 = kind(detector4)
      REAL*8 detector5
      integer, parameter :: cctki_use_detector5 = kind(detector5)
      REAL*8 detector6
      integer, parameter :: cctki_use_detector6 = kind(detector6)
      REAL*8 detector7
      integer, parameter :: cctki_use_detector7 = kind(detector7)
      REAL*8 detector8
      integer, parameter :: cctki_use_detector8 = kind(detector8)
      REAL*8 detector9
      integer, parameter :: cctki_use_detector9 = kind(detector9)
      REAL*8 origin_x
      integer, parameter :: cctki_use_origin_x = kind(origin_x)
      REAL*8 origin_y
      integer, parameter :: cctki_use_origin_y = kind(origin_y)
      REAL*8 origin_z
      integer, parameter :: cctki_use_origin_z = kind(origin_z)
      integer*8 interpolation_operator
      integer, parameter :: cctki_use_interpolation_operator = kind(inte
     &rpolation_operator)
      integer*8 timecoord
      integer, parameter :: cctki_use_timecoord = kind(timecoord)
      INTEGER*4 Cauchy
      integer, parameter :: cctki_use_Cauchy = kind(Cauchy)
      INTEGER*4 Cauchy_timestep
      integer, parameter :: cctki_use_Cauchy_timestep = kind(Cauchy_time
     &step)
      INTEGER*4 Np
      integer, parameter :: cctki_use_Np = kind(Np)
      INTEGER*4 Nt
      integer, parameter :: cctki_use_Nt = kind(Nt)
      INTEGER*4 all_modes
      integer, parameter :: cctki_use_all_modes = kind(all_modes)
      INTEGER*4 doADMmass
      integer, parameter :: cctki_use_doADMmass = kind(doADMmass)
      INTEGER*4 do_momentum
      integer, parameter :: cctki_use_do_momentum = kind(do_momentum)
      INTEGER*4 do_spin
      integer, parameter :: cctki_use_do_spin = kind(do_spin)
      INTEGER*4 interpolation_order
      integer, parameter :: cctki_use_interpolation_order = kind(interpo
     &lation_order)
      INTEGER*4 itout
      integer, parameter :: cctki_use_itout = kind(itout)
      INTEGER*4 l_mode
      integer, parameter :: cctki_use_l_mode = kind(l_mode)
      INTEGER*4 m_mode
      integer, parameter :: cctki_use_m_mode = kind(m_mode)
      INTEGER*4 num_detectors
      integer, parameter :: cctki_use_num_detectors = kind(num_detectors
     &)
      INTEGER*4 verbose
      integer, parameter :: cctki_use_verbose = kind(verbose)
      COMMON /Extractpriv/Cauchy_dr, Cauchy_r1, detector1, detector2, de
     &tector3, detector4, detector5, detector6, detector7, detector8, de
     &tector9, origin_x, origin_y, origin_z, interpolation_operator, tim
     &ecoord, Cauchy, Cauchy_timestep, Np, Nt, all_modes, doADMmass, do_
     &momentum, do_spin, interpolation_order, itout, l_mode, m_mode, num
     &_detectors, verbose
      integer*8 CCTKH0
      integer, parameter :: cctki_use_CCTKH0 = kind(CCTKH0)
      integer*8 CCTKH1
      integer, parameter :: cctki_use_CCTKH1 = kind(CCTKH1)
      integer*8 CCTKH2
      integer, parameter :: cctki_use_CCTKH2 = kind(CCTKH2)
      integer*8 CCTKH3
      integer, parameter :: cctki_use_CCTKH3 = kind(CCTKH3)
      integer*8 CCTKH4
      integer, parameter :: cctki_use_CCTKH4 = kind(CCTKH4)
      integer*8 CCTKH5
      integer, parameter :: cctki_use_CCTKH5 = kind(CCTKH5)
      integer*8 CCTKH6
      integer, parameter :: cctki_use_CCTKH6 = kind(CCTKH6)
      integer*8 CCTKH7
      integer, parameter :: cctki_use_CCTKH7 = kind(CCTKH7)
      integer*8 CCTKH8
      integer, parameter :: cctki_use_CCTKH8 = kind(CCTKH8)
      integer*8 CCTKH9
      integer, parameter :: cctki_use_CCTKH9 = kind(CCTKH9)
      integer*8 CCTKH10
      integer, parameter :: cctki_use_CCTKH10 = kind(CCTKH10)
      integer*8 CCTKH12
      integer, parameter :: cctki_use_CCTKH12 = kind(CCTKH12)
      integer*8 metric_type
      integer, parameter :: cctki_use_metric_type = kind(metric_type)
      integer*8 CCTKH14
      integer, parameter :: cctki_use_CCTKH14 = kind(CCTKH14)
      integer*8 CCTKH15
      integer, parameter :: cctki_use_CCTKH15 = kind(CCTKH15)
      INTEGER*4 CCTKH11
      integer, parameter :: cctki_use_CCTKH11 = kind(CCTKH11)
      INTEGER*4 CCTKH13
      integer, parameter :: cctki_use_CCTKH13 = kind(CCTKH13)
      INTEGER*4 CCTKH16
      integer, parameter :: cctki_use_CCTKH16 = kind(CCTKH16)
      COMMON /ADMBASErest/CCTKH0, CCTKH1, CCTKH2, CCTKH3, CCTKH4, CCTKH5
     &, CCTKH6, CCTKH7, CCTKH8, CCTKH9, CCTKH10, CCTKH12, metric_type, C
     &CTKH14, CCTKH15, CCTKH11, CCTKH13, CCTKH16
      REAL*8 CCTKH18
      integer, parameter :: cctki_use_CCTKH18 = kind(CCTKH18)
      REAL*8 CCTKH19
      integer, parameter :: cctki_use_CCTKH19 = kind(CCTKH19)
      REAL*8 CCTKH20
      integer, parameter :: cctki_use_CCTKH20 = kind(CCTKH20)
      REAL*8 CCTKH21
      integer, parameter :: cctki_use_CCTKH21 = kind(CCTKH21)
      REAL*8 CCTKH31
      integer, parameter :: cctki_use_CCTKH31 = kind(CCTKH31)
      REAL*8 CCTKH32
      integer, parameter :: cctki_use_CCTKH32 = kind(CCTKH32)
      REAL*8 CCTKH33
      integer, parameter :: cctki_use_CCTKH33 = kind(CCTKH33)
      REAL*8 CCTKH34
      integer, parameter :: cctki_use_CCTKH34 = kind(CCTKH34)
      REAL*8 CCTKH35
      integer, parameter :: cctki_use_CCTKH35 = kind(CCTKH35)
      REAL*8 CCTKH36
      integer, parameter :: cctki_use_CCTKH36 = kind(CCTKH36)
      REAL*8 CCTKH37
      integer, parameter :: cctki_use_CCTKH37 = kind(CCTKH37)
      REAL*8 CCTKH38
      integer, parameter :: cctki_use_CCTKH38 = kind(CCTKH38)
      integer*8 CCTKH17
      integer, parameter :: cctki_use_CCTKH17 = kind(CCTKH17)
      integer*8 domain
      integer, parameter :: cctki_use_domain = kind(domain)
      integer*8 CCTKH22
      integer, parameter :: cctki_use_CCTKH22 = kind(CCTKH22)
      integer*8 CCTKH23
      integer, parameter :: cctki_use_CCTKH23 = kind(CCTKH23)
      integer*8 CCTKH30
      integer, parameter :: cctki_use_CCTKH30 = kind(CCTKH30)
      INTEGER*4 CCTKH24
      integer, parameter :: cctki_use_CCTKH24 = kind(CCTKH24)
      INTEGER*4 CCTKH25
      integer, parameter :: cctki_use_CCTKH25 = kind(CCTKH25)
      INTEGER*4 CCTKH26
      integer, parameter :: cctki_use_CCTKH26 = kind(CCTKH26)
      INTEGER*4 CCTKH27
      integer, parameter :: cctki_use_CCTKH27 = kind(CCTKH27)
      INTEGER*4 CCTKH28
      integer, parameter :: cctki_use_CCTKH28 = kind(CCTKH28)
      INTEGER*4 CCTKH29
      integer, parameter :: cctki_use_CCTKH29 = kind(CCTKH29)
      COMMON /GRIDrest/CCTKH18, CCTKH19, CCTKH20, CCTKH21, CCTKH31, CCTK
     &H32, CCTKH33, CCTKH34, CCTKH35, CCTKH36, CCTKH37, CCTKH38, CCTKH17
     &, domain, CCTKH22, CCTKH23, CCTKH30, CCTKH24, CCTKH25, CCTKH26, CC
     &TKH27, CCTKH28, CCTKH29
      REAL*8 CCTKH44
      integer, parameter :: cctki_use_CCTKH44 = kind(CCTKH44)
      REAL*8 CCTKH58
      integer, parameter :: cctki_use_CCTKH58 = kind(CCTKH58)
      REAL*8 CCTKH68
      integer, parameter :: cctki_use_CCTKH68 = kind(CCTKH68)
      REAL*8 CCTKH70
      integer, parameter :: cctki_use_CCTKH70 = kind(CCTKH70)
      REAL*8 CCTKH72
      integer, parameter :: cctki_use_CCTKH72 = kind(CCTKH72)
      REAL*8 CCTKH74
      integer, parameter :: cctki_use_CCTKH74 = kind(CCTKH74)
      REAL*8 CCTKH76
      integer, parameter :: cctki_use_CCTKH76 = kind(CCTKH76)
      REAL*8 CCTKH78
      integer, parameter :: cctki_use_CCTKH78 = kind(CCTKH78)
      REAL*8 CCTKH80
      integer, parameter :: cctki_use_CCTKH80 = kind(CCTKH80)
      REAL*8 CCTKH82
      integer, parameter :: cctki_use_CCTKH82 = kind(CCTKH82)
      REAL*8 CCTKH84
      integer, parameter :: cctki_use_CCTKH84 = kind(CCTKH84)
      integer*8 CCTKH41
      integer, parameter :: cctki_use_CCTKH41 = kind(CCTKH41)
      integer*8 CCTKH42
      integer, parameter :: cctki_use_CCTKH42 = kind(CCTKH42)
      integer*8 CCTKH45
      integer, parameter :: cctki_use_CCTKH45 = kind(CCTKH45)
      integer*8 CCTKH48
      integer, parameter :: cctki_use_CCTKH48 = kind(CCTKH48)
      integer*8 CCTKH49
      integer, parameter :: cctki_use_CCTKH49 = kind(CCTKH49)
      integer*8 CCTKH50
      integer, parameter :: cctki_use_CCTKH50 = kind(CCTKH50)
      integer*8 CCTKH54
      integer, parameter :: cctki_use_CCTKH54 = kind(CCTKH54)
      integer*8 out_dir
      integer, parameter :: cctki_use_out_dir = kind(out_dir)
      integer*8 CCTKH60
      integer, parameter :: cctki_use_CCTKH60 = kind(CCTKH60)
      integer*8 CCTKH61
      integer, parameter :: cctki_use_CCTKH61 = kind(CCTKH61)
      integer*8 CCTKH62
      integer, parameter :: cctki_use_CCTKH62 = kind(CCTKH62)
      integer*8 CCTKH64
      integer, parameter :: cctki_use_CCTKH64 = kind(CCTKH64)
      integer*8 CCTKH86
      integer, parameter :: cctki_use_CCTKH86 = kind(CCTKH86)
      integer*8 CCTKH88
      integer, parameter :: cctki_use_CCTKH88 = kind(CCTKH88)
      integer*8 CCTKH90
      integer, parameter :: cctki_use_CCTKH90 = kind(CCTKH90)
      integer*8 CCTKH92
      integer, parameter :: cctki_use_CCTKH92 = kind(CCTKH92)
      integer*8 CCTKH93
      integer, parameter :: cctki_use_CCTKH93 = kind(CCTKH93)
      integer*8 CCTKH96
      integer, parameter :: cctki_use_CCTKH96 = kind(CCTKH96)
      INTEGER*4 CCTKH39
      integer, parameter :: cctki_use_CCTKH39 = kind(CCTKH39)
      INTEGER*4 CCTKH40
      integer, parameter :: cctki_use_CCTKH40 = kind(CCTKH40)
      INTEGER*4 CCTKH43
      integer, parameter :: cctki_use_CCTKH43 = kind(CCTKH43)
      INTEGER*4 CCTKH46
      integer, parameter :: cctki_use_CCTKH46 = kind(CCTKH46)
      INTEGER*4 CCTKH47
      integer, parameter :: cctki_use_CCTKH47 = kind(CCTKH47)
      INTEGER*4 CCTKH51
      integer, parameter :: cctki_use_CCTKH51 = kind(CCTKH51)
      INTEGER*4 CCTKH52
      integer, parameter :: cctki_use_CCTKH52 = kind(CCTKH52)
      INTEGER*4 CCTKH53
      integer, parameter :: cctki_use_CCTKH53 = kind(CCTKH53)
      INTEGER*4 CCTKH55
      integer, parameter :: cctki_use_CCTKH55 = kind(CCTKH55)
      INTEGER*4 CCTKH56
      integer, parameter :: cctki_use_CCTKH56 = kind(CCTKH56)
      INTEGER*4 CCTKH57
      integer, parameter :: cctki_use_CCTKH57 = kind(CCTKH57)
      INTEGER*4 CCTKH59
      integer, parameter :: cctki_use_CCTKH59 = kind(CCTKH59)
      INTEGER*4 CCTKH63
      integer, parameter :: cctki_use_CCTKH63 = kind(CCTKH63)
      INTEGER*4 CCTKH65
      integer, parameter :: cctki_use_CCTKH65 = kind(CCTKH65)
      INTEGER*4 CCTKH66
      integer, parameter :: cctki_use_CCTKH66 = kind(CCTKH66)
      INTEGER*4 CCTKH67
      integer, parameter :: cctki_use_CCTKH67 = kind(CCTKH67)
      INTEGER*4 CCTKH69
      integer, parameter :: cctki_use_CCTKH69 = kind(CCTKH69)
      INTEGER*4 CCTKH71
      integer, parameter :: cctki_use_CCTKH71 = kind(CCTKH71)
      INTEGER*4 CCTKH73
      integer, parameter :: cctki_use_CCTKH73 = kind(CCTKH73)
      INTEGER*4 CCTKH75
      integer, parameter :: cctki_use_CCTKH75 = kind(CCTKH75)
      INTEGER*4 CCTKH77
      integer, parameter :: cctki_use_CCTKH77 = kind(CCTKH77)
      INTEGER*4 CCTKH79
      integer, parameter :: cctki_use_CCTKH79 = kind(CCTKH79)
      INTEGER*4 CCTKH81
      integer, parameter :: cctki_use_CCTKH81 = kind(CCTKH81)
      INTEGER*4 CCTKH83
      integer, parameter :: cctki_use_CCTKH83 = kind(CCTKH83)
      INTEGER*4 CCTKH85
      integer, parameter :: cctki_use_CCTKH85 = kind(CCTKH85)
      INTEGER*4 CCTKH87
      integer, parameter :: cctki_use_CCTKH87 = kind(CCTKH87)
      INTEGER*4 CCTKH89
      integer, parameter :: cctki_use_CCTKH89 = kind(CCTKH89)
      INTEGER*4 CCTKH91
      integer, parameter :: cctki_use_CCTKH91 = kind(CCTKH91)
      INTEGER*4 CCTKH94
      integer, parameter :: cctki_use_CCTKH94 = kind(CCTKH94)
      INTEGER*4 CCTKH95
      integer, parameter :: cctki_use_CCTKH95 = kind(CCTKH95)
      COMMON /IOrest/CCTKH44, CCTKH58, CCTKH68, CCTKH70, CCTKH72, CCTKH7
     &4, CCTKH76, CCTKH78, CCTKH80, CCTKH82, CCTKH84, CCTKH41, CCTKH42, 
     &CCTKH45, CCTKH48, CCTKH49, CCTKH50, CCTKH54, out_dir, CCTKH60, CCT
     &KH61, CCTKH62, CCTKH64, CCTKH86, CCTKH88, CCTKH90, CCTKH92, CCTKH9
     &3, CCTKH96, CCTKH39, CCTKH40, CCTKH43, CCTKH46, CCTKH47, CCTKH51, 
     &CCTKH52, CCTKH53, CCTKH55, CCTKH56, CCTKH57, CCTKH59, CCTKH63, CCT
     &KH65, CCTKH66, CCTKH67, CCTKH69, CCTKH71, CCTKH73, CCTKH75, CCTKH7
     &7, CCTKH79, CCTKH81, CCTKH83, CCTKH85, CCTKH87, CCTKH89, CCTKH91, 
     &CCTKH94, CCTKH95
      interface
      integer*8 function CCTK_PointerTo (var)
      implicit none
      type(*),dimension(..),TARGET :: var
      end function CCTK_PointerTo
      end interface
      interface
      integer function CCTK_Equals (arg1, arg2)
      implicit none
      integer*8 arg1
      character(*) arg2
      end function CCTK_Equals
      integer function CCTK_MyProc (cctkGH)
      implicit none
      integer*8 cctkGH
      end function CCTK_MyProc
      integer function CCTK_nProcs (cctkGH)
      implicit none
      integer*8 cctkGH
      end function CCTK_nProcs
      integer function CCTK_IsThornActive (name)
      implicit none
      character(*) name
      end function CCTK_IsThornActive
      integer*8 function CCTK_NullPointer ()
      implicit none
      end function CCTK_NullPointer
      end interface

c

      Pi = ACOS(-1D0)

c
      SELECT CASE (Psi_power)

      CASE (1)
         ip = 4

      CASE (4)
         ip = 1

      CASE DEFAULT
         WRITE(*,*) "This value of Psi_power is not supported"

      END SELECT   
      

      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

                  ux = (x(i)-origin(1))/rad
                  uy = (y(j)-origin(2))/rad
                  uz = (z(k)-origin(3))/rad

c
c

                  if (conformal_state > 0) then
  
                    p = psi(i,j,k)**ip

                    dxx = p*gxx(i,j,k); dxy = p*gxy(i,j,k)
                    dxz = p*gxz(i,j,k); dyy = p*gyy(i,j,k)
                    dyz = p*gyz(i,j,k); dzz = p*gzz(i,j,k)

                 else 

                    p = 1.0d0

                    dxx = gxx(i,j,k); dxy = gxy(i,j,k)
                    dxz = gxz(i,j,k); dyy = gyy(i,j,k)
                    dyz = gyz(i,j,k); dzz = gzz(i,j,k)

                 end if

c
c
                  det = (dxx*dyy*dzz + 2.0D0*dxy*dxz*dyz
     &                 - (dxx*dyz**2 + dyy*dxz**2 + dzz*dxy**2))

                  idet = 1.0/det

c
c
                  uxx = idet*(dyy*dzz - dyz**2)
                  uyy = idet*(dxx*dzz - dxz**2)
                  uzz = idet*(dxx*dyy - dxy**2)
                  uxy = idet*(dxz*dyz - dzz*dxy)
                  uxz = idet*(dxy*dyz - dyy*dxz)
                  uyz = idet*(dxy*dxz - dxx*dyz)

c
c

                  IF (conformal_state > 0) THEN

                    pip = psi(i+1,j,k)**ip
                    pim = psi(i-1,j,k)**ip
                    pjp = psi(i,j+1,k)**ip
                    pjm = psi(i,j-1,k)**ip
                    pkp = psi(i,j,k+1)**ip
                    pkm = psi(i,j,k-1)**ip

                 else

                    pip = 1.0d0
                    pim = 1.0d0
                    pjp = 1.0d0
                    pjm = 1.0d0
                    pkp = 1.0d0
                    pkm = 1.0d0

                 end if

                  dxx_y = half/Dy*(pjp*gxx(i,j+1,k)-pjm*gxx(i,j-1,k))
                  dxx_z = half/Dz*(pkp*gxx(i,j,k+1)-pkm*gxx(i,j,k-1))
                  dxy_x = half/Dx*(pip*gxy(i+1,j,k)-pim*gxy(i-1,j,k))
                  dxy_y = half/Dy*(pjp*gxy(i,j+1,k)-pjm*gxy(i,j-1,k))
                  dxy_z = half/Dz*(pkp*gxy(i,j,k+1)-pkm*gxy(i,j,k-1))
                  dyy_x = half/Dx*(pip*gyy(i+1,j,k)-pim*gyy(i-1,j,k))
                  dyy_z = half/Dz*(pkp*gyy(i,j,k+1)-pkm*gyy(i,j,k-1))
                  dxz_x = half/Dx*(pip*gxz(i+1,j,k)-pim*gxz(i-1,j,k))
                  dxz_y = half/Dy*(pjp*gxz(i,j+1,k)-pjm*gxz(i,j-1,k))
                  dxz_z = half/Dz*(pkp*gxz(i,j,k+1)-pkm*gxz(i,j,k-1))
                  dyz_x = half/Dx*(pip*gyz(i+1,j,k)-pim*gyz(i-1,j,k))
                  dyz_y = half/Dy*(pjp*gyz(i,j+1,k)-pjm*gyz(i,j-1,k))
                  dyz_z = half/Dz*(pkp*gyz(i,j,k+1)-pkm*gyz(i,j,k-1))
                  dzz_x = half/Dx*(pip*gzz(i+1,j,k)-pim*gzz(i-1,j,k))
                  dzz_y = half/Dy*(pjp*gzz(i,j+1,k)-pjm*gzz(i,j-1,k))

                  term1 = uxy*(dxx_y-dxy_x)+uxz*(dxx_z-dxz_x)
     &                 +uyy*(dxy_y-dyy_x)+uyz*(dxy_z-dyz_x)
     &                 +uyz*(dxz_y-dyz_x)+uzz*(dxz_z-dzz_x)

                  term2 = uyz*(dyy_z-dyz_y)+uxy*(dyy_x-dxy_y)
     &                 +uzz*(dyz_z-dzz_y)+uxz*(dyz_x-dxz_y)
     &                 +uxz*(dxy_z-dxz_y)+uxx*(dxy_x-dxx_y)

                  term3 = uxz*(dzz_x-dxz_z)+uyz*(dzz_y-dyz_z)
     &                 +uxx*(dxz_x-dxx_z)+uxy*(dxz_y-dxy_z)
     &                 +uxy*(dyz_x-dxy_z)+uyy*(dyz_y-dyy_z)

                  ADMmass_int(i,j,k) = 1.0D0/16.0D0/Pi*(ux*term1+
     &                 uy*term2+uz*term3)*SQRT(det)*rad**2

               ELSE
                  
                  ADMmass_int(i,j,k) = 0.0D0
                  
               ENDIF

          ENDDO
        ENDDO
      ENDDO


c
c

      DO k = 2, size(z)-1
         DO j = 2, size(y)-1
            ADMmass_int(1,j,k) = ADMmass_int(2,j,k)
         ENDDO
      ENDDO
      DO k = 2, size(z)-1
         DO i = 1, size(x)-1
            ADMmass_int(i,1,k) = ADMmass_int(i,2,k)
         ENDDO
      ENDDO
      DO j = 1, size(y)-1
         DO i = 1, size(x)-1
            ADMmass_int(i,j,1) = ADMmass_int(i,j,2)
         ENDDO
      ENDDO


      END SUBROUTINE ADMmass_integrand3D
