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

      subroutine EvolSimple_KSources(CCTK_ARGUMENTS)

      implicit none

      DECLARE_CCTK_ARGUMENTS
      DECLARE_CCTK_PARAMETERS
      DECLARE_CCTK_FUNCTIONS

      CCTK_REAL detg_tempxx,detg_tempxy,detg_tempxz 
      CCTK_REAL detg_tempyy,detg_tempyz,detg_tempzz 
      CCTK_REAL detg_detg,detg_detcg 
      CCTK_REAL uppermet_fdet 
      CCTK_REAL uppermet_uxx,uppermet_uxy,uppermet_uxz 
      CCTK_REAL uppermet_uyy,uppermet_uyz,uppermet_uzz 
      CCTK_REAL delgb111,delgb112,delgb113,delgb122,delgb123,delgb133 
      CCTK_REAL dxdgxx,dxdgxy,dxdgxz,dxdgyy,dxdgyz,dxdgzz
      CCTK_REAL dxxdgxx,dxxdgxy,dxxdgxz,dxxdgyy,dxxdgyz,dxxdgzz 
      CCTK_REAL delgb211,delgb212,delgb213,delgb222,delgb223,delgb233 
      CCTK_REAL dydgxx,dydgxy,dydgxz,dydgyy,dydgyz,dydgzz
      CCTK_REAL deldelg1211,deldelg1212,deldelg1213,deldelg1222,deldelg1223,deldelg1233 
      CCTK_REAL delgb311,delgb312,delgb313,delgb322,delgb323,delgb333 
      CCTK_REAL dzdgxx,dzdgxy,dzdgxz,dzdgyy,dzdgyz,dzdgzz
      CCTK_REAL deldelg1311,deldelg1312,deldelg1313,deldelg1322,deldelg1323,deldelg1333 
      CCTK_REAL dyydgxx,dyydgxy,dyydgxz 
      CCTK_REAL dyydgyy,dyydgyz,dyydgzz 
      CCTK_REAL deldelg2311,deldelg2312,deldelg2313,deldelg2322,deldelg2323,deldelg2333 
      CCTK_REAL dzzdgxx,dzzdgxy,dzzdgxz,dzzdgyy,dzzdgyz,dzzdgzz 
      CCTK_REAL gammado111,gammado112,gammado113,gammado122,gammado123,gammado133 
      CCTK_REAL gammado211,gammado212,gammado213,gammado222,gammado223,gammado233 
      CCTK_REAL gammado311,gammado312,gammado313,gammado322,gammado323,gammado333 
      CCTK_REAL gamma111,gamma112,gamma113,gamma122 
      CCTK_REAL gamma123,gamma133,gamma211,gamma212,gamma213,gamma222,gamma223 
      CCTK_REAL gamma233,gamma311,gamma312,gamma313,gamma322,gamma323,gamma333 
      CCTK_REAL ricci_R11,ricci_R12,ricci_R13,ricci_R22,ricci_R23,ricci_R33
      CCTK_REAL KK11,KK12,KK13,KK22,KK23,KK33
      CCTK_REAL trk_trK 
      CCTK_REAL da_dxda,da_dyda,da_dzda 
      CCTK_REAL oodx2,oody2,oodz2
      CCTK_REAL oo2dx,oo2dy,oo2dz 
      CCTK_REAL oo4dxdy,oo4dxdz,oo4dydz 
      CCTK_REAL dda_dxxda,dda_dxyda,dda_dxzda,dda_dyyda,dda_dyzda,dda_dzzda 
      CCTK_REAL cdcda_cdxxda,cdcda_cdxyda,cdcda_cdxzda,cdcda_cdyyda,cdcda_cdyzda,cdcda_cdzzda

      integer i,j,k
      integer nx,ny,nz

      CCTK_REAL :: dx,dy,dz,dt

c     Grid.

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

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

      oo4dydz = 1D0/(4D0*dy*dz)
      oo4dxdy = 1D0/(4D0*dx*dy)
      oo4dxdz = 1D0/(4D0*dx*dz)

      oo2dx = 1D0/(2D0*dx)
      oo2dy = 1D0/(2D0*dy)
      oo2dz = 1D0/(2D0*dz)

      oodx2 = 1D0/(dx*dx)
      oody2 = 1D0/(dy*dy)
      oodz2 = 1D0/(dz*dz)

      do k=2,nz-1
         do j=2,ny-1
            do i=2,nx-1

               detg_tempxx =  gyy_p(i,j,k)*gzz_p(i,j,k)-gyz_p(i,j,k)*gyz_p(i,j,k)
               detg_tempxy =  gxz_p(i,j,k)*gyz_p(i,j,k)-gxy_p(i,j,k)*gzz_p(i,j,k)
               detg_tempxz = -gxz_p(i,j,k)*gyy_p(i,j,k)+gxy_p(i,j,k)*gyz_p(i,j,k) 
               detg_tempyy =  gxx_p(i,j,k)*gzz_p(i,j,k)-gxz_p(i,j,k)*gxz_p(i,j,k)
               detg_tempyz =  gxy_p(i,j,k)*gxz_p(i,j,k)-gxx_p(i,j,k)*gyz_p(i,j,k)
               detg_tempzz =  gxx_p(i,j,k)*gyy_p(i,j,k)-gxy_p(i,j,k)*gxy_p(i,j,k)

               detg_detcg = (detg_tempxx*gxx_p(i,j,k)+detg_tempxy*gxy_p(i,j,k)+detg_tempxz*gxz_p(i,j,k))

               uppermet_fdet = detg_detcg 

               uppermet_uxx = detg_tempxx / uppermet_fdet 
               uppermet_uxy = detg_tempxy / uppermet_fdet 
               uppermet_uxz = detg_tempxz / uppermet_fdet 
               uppermet_uyy = detg_tempyy / uppermet_fdet 
               uppermet_uyz = detg_tempyz / uppermet_fdet 
               uppermet_uzz = detg_tempzz / uppermet_fdet 

               delgb111 = oo2dx*(gxx_p(i+1,j,k)-gxx_p(i-1,j,k))
               delgb112 = oo2dx*(gxy_p(i+1,j,k)-gxy_p(i-1,j,k))
               delgb113 = oo2dx*(gxz_p(i+1,j,k)-gxz_p(i-1,j,k))
               delgb122 = oo2dx*(gyy_p(i+1,j,k)-gyy_p(i-1,j,k))
               delgb123 = oo2dx*(gyz_p(i+1,j,k)-gyz_p(i-1,j,k))
               delgb133 = oo2dx*(gzz_p(i+1,j,k)-gzz_p(i-1,j,k))

               dxdgxx = delgb111
               dxdgxy = delgb112
               dxdgxz = delgb113
               dxdgyy = delgb122
               dxdgyz = delgb123
               dxdgzz = delgb133

               dxxdgxx = oodx2*(gxx_p(i+1,j,k)-2*gxx_p(i,j,k)+gxx_p(i-1,j,k))
               dxxdgxy = oodx2*(gxy_p(i+1,j,k)-2*gxy_p(i,j,k)+gxy_p(i-1,j,k))
               dxxdgxz = oodx2*(gxz_p(i+1,j,k)-2*gxz_p(i,j,k)+gxz_p(i-1,j,k))
               dxxdgyy = oodx2*(gyy_p(i+1,j,k)-2*gyy_p(i,j,k)+gyy_p(i-1,j,k))
               dxxdgyz = oodx2*(gyz_p(i+1,j,k)-2*gyz_p(i,j,k)+gyz_p(i-1,j,k))
               dxxdgzz = oodx2*(gzz_p(i+1,j,k)-2*gzz_p(i,j,k)+gzz_p(i-1,j,k))

               delgb211 = oo2dy*(gxx_p(i,j+1,k)-gxx_p(i,j-1,k))
               delgb212 = oo2dy*(gxy_p(i,j+1,k)-gxy_p(i,j-1,k))
               delgb213 = oo2dy*(gxz_p(i,j+1,k)-gxz_p(i,j-1,k))
               delgb222 = oo2dy*(gyy_p(i,j+1,k)-gyy_p(i,j-1,k))
               delgb223 = oo2dy*(gyz_p(i,j+1,k)-gyz_p(i,j-1,k))
               delgb233 = oo2dy*(gzz_p(i,j+1,k)-gzz_p(i,j-1,k))

               dydgxx = delgb211
               dydgxy = delgb212
               dydgxz = delgb213
               dydgyy = delgb222
               dydgyz = delgb223
               dydgzz = delgb233

               deldelg1211 = oo4dxdy*(gxx_p(i+1,j+1,k)-gxx_p(i+1,j-1,k)-gxx_p(i-1,j+1,k)+gxx_p(i-1,j-1,k))
               deldelg1212 = oo4dxdy*(gxy_p(i+1,j+1,k)-gxy_p(i+1,j-1,k)-gxy_p(i-1,j+1,k)+gxy_p(i-1,j-1,k))
               deldelg1213 = oo4dxdy*(gxz_p(i+1,j+1,k)-gxz_p(i+1,j-1,k)-gxz_p(i-1,j+1,k)+gxz_p(i-1,j-1,k))
               deldelg1222 = oo4dxdy*(gyy_p(i+1,j+1,k)-gyy_p(i+1,j-1,k)-gyy_p(i-1,j+1,k)+gyy_p(i-1,j-1,k))
               deldelg1223 = oo4dxdy*(gyz_p(i+1,j+1,k)-gyz_p(i+1,j-1,k)-gyz_p(i-1,j+1,k)+gyz_p(i-1,j-1,k))
               deldelg1233 = oo4dxdy*(gzz_p(i+1,j+1,k)-gzz_p(i+1,j-1,k)-gzz_p(i-1,j+1,k)+gzz_p(i-1,j-1,k))

               delgb311 = oo2dz*(gxx_p(i,j,k+1)-gxx_p(i,j,k-1))
               delgb312 = oo2dz*(gxy_p(i,j,k+1)-gxy_p(i,j,k-1))
               delgb313 = oo2dz*(gxz_p(i,j,k+1)-gxz_p(i,j,k-1))
               delgb322 = oo2dz*(gyy_p(i,j,k+1)-gyy_p(i,j,k-1))
               delgb323 = oo2dz*(gyz_p(i,j,k+1)-gyz_p(i,j,k-1))
               delgb333 = oo2dz*(gzz_p(i,j,k+1)-gzz_p(i,j,k-1))

               dzdgxx = delgb311
               dzdgxy = delgb312
               dzdgxz = delgb313
               dzdgyy = delgb322
               dzdgyz = delgb323
               dzdgzz = delgb333

               deldelg1311 = oo4dxdz*(gxx_p(i+1,j,k+1)-gxx_p(i+1,j,k-1)-gxx_p(i-1,j,k+1)+gxx_p(i-1,j,k-1))
               deldelg1312 = oo4dxdz*(gxy_p(i+1,j,k+1)-gxy_p(i+1,j,k-1)-gxy_p(i-1,j,k+1)+gxy_p(i-1,j,k-1))
               deldelg1313 = oo4dxdz*(gxz_p(i+1,j,k+1)-gxz_p(i+1,j,k-1)-gxz_p(i-1,j,k+1)+gxz_p(i-1,j,k-1))
               deldelg1322 = oo4dxdz*(gyy_p(i+1,j,k+1)-gyy_p(i+1,j,k-1)-gyy_p(i-1,j,k+1)+gyy_p(i-1,j,k-1))
               deldelg1323 = oo4dxdz*(gyz_p(i+1,j,k+1)-gyz_p(i+1,j,k-1)-gyz_p(i-1,j,k+1)+gyz_p(i-1,j,k-1))
               deldelg1333 = oo4dxdz*(gzz_p(i+1,j,k+1)-gzz_p(i+1,j,k-1)-gzz_p(i-1,j,k+1)+gzz_p(i-1,j,k-1))

               dyydgxx = oody2*(gxx_p(i,j+1,k)-2*gxx_p(i,j,k)+gxx_p(i,j-1,k))
               dyydgxy = oody2*(gxy_p(i,j+1,k)-2*gxy_p(i,j,k)+gxy_p(i,j-1,k))
               dyydgxz = oody2*(gxz_p(i,j+1,k)-2*gxz_p(i,j,k)+gxz_p(i,j-1,k))
               dyydgyy = oody2*(gyy_p(i,j+1,k)-2*gyy_p(i,j,k)+gyy_p(i,j-1,k))
               dyydgyz = oody2*(gyz_p(i,j+1,k)-2*gyz_p(i,j,k)+gyz_p(i,j-1,k))
               dyydgzz = oody2*(gzz_p(i,j+1,k)-2*gzz_p(i,j,k)+gzz_p(i,j-1,k))
 
               deldelg2311 = oo4dydz*(gxx_p(i,j+1,k+1)-gxx_p(i,j+1,k-1)-gxx_p(i,j-1,k+1)+gxx_p(i,j-1,k-1))
               deldelg2312 = oo4dydz*(gxy_p(i,j+1,k+1)-gxy_p(i,j+1,k-1)-gxy_p(i,j-1,k+1) +gxy_p(i,j-1,k-1))
               deldelg2313 = oo4dydz*(gxz_p(i,j+1,k+1)-gxz_p(i,j+1,k-1)-gxz_p(i,j-1,k+1)+gxz_p(i,j-1,k-1))
               deldelg2322 = oo4dydz*(gyy_p(i,j+1,k+1)-gyy_p(i,j+1,k-1)-gyy_p(i,j-1,k+1)+gyy_p(i,j-1,k-1))
               deldelg2323 = oo4dydz*(gyz_p(i,j+1,k+1)-gyz_p(i,j+1,k-1)-gyz_p(i,j-1,k+1)+gyz_p(i,j-1,k-1))
               deldelg2333 = oo4dydz*(gzz_p(i,j+1,k+1)-gzz_p(i,j+1,k-1)-gzz_p(i,j-1,k+1)+gzz_p(i,j-1,k-1))

               dzzdgxx = oodz2*(gxx_p(i,j,k+1)-2*gxx_p(i,j,k)+gxx_p(i,j,k-1))
               dzzdgxy = oodz2*(gxy_p(i,j,k+1)-2*gxy_p(i,j,k)+gxy_p(i,j,k-1))
               dzzdgxz = oodz2*(gxz_p(i,j,k+1)-2*gxz_p(i,j,k)+gxz_p(i,j,k-1))
               dzzdgyy = oodz2*(gyy_p(i,j,k+1)-2*gyy_p(i,j,k)+gyy_p(i,j,k-1))
               dzzdgyz = oodz2*(gyz_p(i,j,k+1)-2*gyz_p(i,j,k)+gyz_p(i,j,k-1))
               dzzdgzz = oodz2*(gzz_p(i,j,k+1)-2*gzz_p(i,j,k)+gzz_p(i,j,k-1))

               gammado111 = dxdgxx/2D0
               gammado112 = dydgxx/2D0
               gammado113 = dzdgxx/2D0
               gammado122 = -dxdgyy/2D0+dydgxy 
               gammado123 = (-dxdgyz+dydgxz+dzdgxy)/2D0
               gammado133 = -dxdgzz/2D0+dzdgxz 
               gammado211 = dxdgxy-dydgxx/2D0
               gammado212 = dxdgyy/2D0
               gammado213 = (dxdgyz-dydgxz+dzdgxy)/2D0
               gammado222 = dydgyy/2D0
               gammado223 = dzdgyy/2D0
               gammado233 = -dydgzz/2D0+dzdgyz 
               gammado311 = dxdgxz-dzdgxx/2D0
               gammado312 = (dxdgyz+dydgxz-dzdgxy)/2D0
               gammado313 = dxdgzz/2D0
               gammado322 = dydgyz-dzdgyy/2D0
               gammado323 = dydgzz/2D0
               gammado333 = dzdgzz/2D0

               gamma111 = gammado111*uppermet_uxx+gammado211*uppermet_uxy+gammado311*uppermet_uxz 
               gamma112 = gammado112*uppermet_uxx+gammado212*uppermet_uxy+gammado312*uppermet_uxz 
               gamma113 = gammado113*uppermet_uxx+gammado213*uppermet_uxy+gammado313*uppermet_uxz 
               gamma122 = gammado122*uppermet_uxx+gammado222*uppermet_uxy+gammado322*uppermet_uxz 
               gamma123 = gammado123*uppermet_uxx+gammado223*uppermet_uxy+gammado323*uppermet_uxz 
               gamma133 = gammado133*uppermet_uxx+gammado233*uppermet_uxy+gammado333*uppermet_uxz 
               gamma211 = gammado111*uppermet_uxy+gammado211*uppermet_uyy+gammado311*uppermet_uyz 
               gamma212 = gammado112*uppermet_uxy+gammado212*uppermet_uyy+gammado312*uppermet_uyz 
               gamma213 = gammado113*uppermet_uxy+gammado213*uppermet_uyy+gammado313*uppermet_uyz 
               gamma222 = gammado122*uppermet_uxy+gammado222*uppermet_uyy+gammado322*uppermet_uyz 
               gamma223 = gammado123*uppermet_uxy+gammado223*uppermet_uyy+gammado323*uppermet_uyz 
               gamma233 = gammado133*uppermet_uxy+gammado233*uppermet_uyy+gammado333*uppermet_uyz 
               gamma311 = gammado111*uppermet_uxz+gammado211*uppermet_uyz+gammado311*uppermet_uzz 
               gamma312 = gammado112*uppermet_uxz+gammado212*uppermet_uyz+gammado312*uppermet_uzz 
               gamma313 = gammado113*uppermet_uxz+gammado213*uppermet_uyz+gammado313*uppermet_uzz 
               gamma322 = gammado122*uppermet_uxz+gammado222*uppermet_uyz+gammado322*uppermet_uzz 
               gamma323 = gammado123*uppermet_uxz+gammado223*uppermet_uyz+gammado323*uppermet_uzz 
               gamma333 = gammado133*uppermet_uxz+gammado233*uppermet_uyz+gammado333*uppermet_uzz 

               ricci_R11 = (deldelg1212+(-dxxdgyy-dyydgxx)/2-
     &              gammado122*gamma111+gammado112*gamma112-gammado222*gamma211 +
     &              gammado212*gamma212-gammado322*gamma311 +
     &              gammado312*gamma312)*uppermet_uyy+(-dxxdgyz+deldelg1213 +
     &              deldelg1312-deldelg2311-2*(gammado223*gamma211 +
     &              gammado323*gamma311)+2*(-(gammado123*gamma111)+gammado113*gamma112 +
     &              gammado213*gamma212+gammado313*gamma312))*uppermet_uyz+(deldelg1313 
     &              +(-dxxdgzz-dzzdgxx)/2-gammado133*gamma111 +
     &              gammado113*gamma113-gammado233*gamma211+gammado213*gamma213 -
     &              gammado333*gamma311+gammado313*gamma313)*uppermet_uzz 
               ricci_R12 = (-deldelg1212+(dxxdgyy+dyydgxx)/2+
     &              gammado122*gamma111-gammado112*gamma112+gammado222*gamma211 -
     &              gammado212*gamma212+gammado322*gamma311 -
     &              gammado312*gamma312)*uppermet_uxy+((dxxdgyz-deldelg1213-
     &              deldelg1312+deldelg2311)/2+gammado123*gamma111 -
     &              gammado113*gamma112+gammado223*gamma211-gammado213*gamma212 +
     &              gammado323*gamma311-gammado313*gamma312)*uppermet_uxz +
     &              ((-deldelg1223+deldelg1322+dyydgxz-deldelg2312)/2-
     &              gammado123*gamma112+gammado122*gamma113-gammado223*gamma212 +
     &              gammado222*gamma213-gammado323*gamma312 +
     &              gammado322*gamma313)*uppermet_uyz+((-deldelg1233+deldelg1323 +
     &              deldelg2313-dzzdgxy)/2-gammado133*gamma112 +
     &              gammado123*gamma113-gammado233*gamma212+gammado223*gamma213 -
     &              gammado333*gamma312+gammado323*gamma313)*uppermet_uzz
               ricci_R13 = ((dxxdgyz-deldelg1213-deldelg1312 +
     &              deldelg2311)/2+gammado123*gamma111-gammado113*gamma112+
     &              gammado223*gamma211-gammado213*gamma212+gammado323*gamma311 -
     &              gammado313*gamma312)*uppermet_uxy+(-deldelg1313+(dxxdgzz+
     &              dzzdgxx)/2+gammado133*gamma111-gammado113*gamma113 +
     &              gammado233*gamma211-gammado213*gamma213+gammado333*gamma311 -
     &              gammado313*gamma313)*uppermet_uxz+((deldelg1223-deldelg1322 -
     &              dyydgxz+deldelg2312)/2+gammado123*gamma112 -
     &              gammado122*gamma113+gammado223*gamma212-gammado222*gamma213 +
     &              gammado323*gamma312-gammado322*gamma313)*uppermet_uyy+((deldelg1233 
     &              -deldelg1323-deldelg2313+dzzdgxy)/2+gammado133*gamma112 
     &              -gammado123*gamma113+gammado233*gamma212-gammado223*gamma213+
     &              gammado333*gamma312-gammado323*gamma313)*uppermet_uyz
               ricci_R22 = (deldelg1212+(-dxxdgyy-dyydgxx)/2-
     &              gammado122*gamma111+gammado112*gamma112-gammado222*gamma211 +
     &              gammado212*gamma212-gammado322*gamma311 +
     &              gammado312*gamma312)*uppermet_uxx+(deldelg1223-deldelg1322 -
     &              dyydgxz+deldelg2312+2*(gammado123*gamma112 -
     &              gammado122*gamma113+gammado223*gamma212+gammado323*gamma312)-
     &              2*(gammado222*gamma213+gammado322*gamma313))*uppermet_uxz +
     &              (deldelg2323+(-dyydgzz-dzzdgyy)/2-
     &              gammado133*gamma122+gammado123*gamma123-gammado233*gamma222 +
     &              gammado223*gamma223-gammado333*gamma322 +
     &              gammado323*gamma323)*uppermet_uzz
               ricci_R23 = ((-dxxdgyz+deldelg1213+deldelg1312 -
     &              deldelg2311)/2-gammado123*gamma111+gammado113*gamma112-
     &              gammado223*gamma211+gammado213*gamma212-gammado323*gamma311 +
     &              gammado313*gamma312)*uppermet_uxx+((-deldelg1223+deldelg1322 +
     &              dyydgxz-deldelg2312)/2-gammado123*gamma112 +
     &              gammado122*gamma113-gammado223*gamma212+gammado222*gamma213 -
     &              gammado323*gamma312+gammado322*gamma313)*uppermet_uxy+((deldelg1233 
     &              -deldelg1323-deldelg2313+dzzdgxy)/2+gammado133*gamma112 
     &              -gammado123*gamma113+gammado233*gamma212-gammado223*gamma213 +
     &              gammado333*gamma312-gammado323*gamma313)*uppermet_uxz+(-deldelg2323 
     &              +(dyydgzz+dzzdgyy)/2+gammado133*gamma122 -
     &              gammado123*gamma123+gammado233*gamma222-gammado223*gamma223 +
     &              gammado333*gamma322-gammado323*gamma323)*uppermet_uyz
               ricci_R33 = (deldelg1313+(-dxxdgzz-dzzdgxx)/2-
     &              gammado133*gamma111+gammado113*gamma113-gammado233*gamma211 +
     &              gammado213*gamma213-gammado333*gamma311 +
     &              gammado313*gamma313)*uppermet_uxx+(-deldelg1233+deldelg1323 +
     &              deldelg2313-dzzdgxy-2*(gammado233*gamma212 +
     &              gammado333*gamma312)+2*(-(gammado133*gamma112)+gammado123*gamma113 +
     &              gammado223*gamma213+gammado323*gamma313))*uppermet_uxy  
     &              + (deldelg2323+(-dyydgzz-dzzdgyy)/2-gammado133*gamma122 +
     &              gammado123*gamma123-gammado233*gamma222+gammado223*gamma223 -
     &              gammado333*gamma322+gammado323*gamma323)*uppermet_uyy

               KK11 = 2*(uppermet_uyz*kxy_p(i,j,k)*kxz_p(i,j,k)+kxx_p(i,j,k)*(uppermet_uxy*kxy_p(i,j,k) +
     &              uppermet_uxz*kxz_p(i,j,k)))+uppermet_uxx*kxx_p(i,j,k)*kxx_p(i,j,k) +
     &              uppermet_uyy*kxy_p(i,j,k)*kxy_p(i,j,k)+uppermet_uzz*kxz_p(i,j,k)*kxz_p(i,j,k) 
               KK12 = kxy_p(i,j,k)*(uppermet_uxx*kxx_p(i,j,k)+uppermet_uyy*kyy_p(i,j,k))+
     &              uppermet_uzz*kxz_p(i,j,k)*kyz_p(i,j,k)+uppermet_uxz*(kxy_p(i,j,k)*kxz_p(i,j,k) +
     &              kxx_p(i,j,k)*kyz_p(i,j,k))+uppermet_uyz*(kxz_p(i,j,k)*kyy_p(i,j,k)+kxy_p(i,j,k)*kyz_p(i,j,k))+
     &              uppermet_uxy*(kxx_p(i,j,k)*kyy_p(i,j,k)+kxy_p(i,j,k)*kxy_p(i,j,k)) 
               KK13 = uppermet_uyy*kxy_p(i,j,k)*kyz_p(i,j,k)+uppermet_uxy*(kxy_p(i,j,k)*kxz_p(i,j,k) +
     &              kxx_p(i,j,k)*kyz_p(i,j,k))+kxz_p(i,j,k)*(uppermet_uxx*kxx_p(i,j,k)+uppermet_uzz*kzz_p(i,j,k))+
     &              uppermet_uyz*(kxz_p(i,j,k)*kyz_p(i,j,k)+kxy_p(i,j,k)*kzz_p(i,j,k))+
     &              uppermet_uxz*(kxx_p(i,j,k)*kzz_p(i,j,k)+kxz_p(i,j,k)*kxz_p(i,j,k)) 
               KK22 = 2*(uppermet_uyz*kyy_p(i,j,k)*kyz_p(i,j,k)+kxy_p(i,j,k)*(uppermet_uxy*kyy_p(i,j,k) +
     &              uppermet_uxz*kyz_p(i,j,k)))+uppermet_uxx*kxy_p(i,j,k)*kxy_p(i,j,k)  +
     &              uppermet_uyy*kyy_p(i,j,k)*kyy_p(i,j,k)+uppermet_uzz*kyz_p(i,j,k)*kyz_p(i,j,k)
               KK23 = uppermet_uxx*kxy_p(i,j,k)*kxz_p(i,j,k)+uppermet_uxy*(kxz_p(i,j,k)*kyy_p(i,j,k) +
     &              kxy_p(i,j,k)*kyz_p(i,j,k))+kyz_p(i,j,k)*(uppermet_uyy*kyy_p(i,j,k)+uppermet_uzz*kzz_p(i,j,k))+
     &              uppermet_uxz*(kxz_p(i,j,k)*kyz_p(i,j,k)+kxy_p(i,j,k)*kzz_p(i,j,k))+
     &              uppermet_uyz*(kyy_p(i,j,k)*kzz_p(i,j,k)+kyz_p(i,j,k)*kyz_p(i,j,k)) 
               KK33 = 2*(uppermet_uyz*kyz_p(i,j,k)*kzz_p(i,j,k)+kxz_p(i,j,k)*(
     &              uppermet_uxy*kyz_p(i,j,k) +
     &              uppermet_uxz*kzz_p(i,j,k)))+uppermet_uxx*kxz_p(i,j,k)*kxz_p(i,j,k) +
     &              uppermet_uyy*kyz_p(i,j,k)*kyz_p(i,j,k)+uppermet_uzz*kzz_p(i,j,k)*kzz_p(i,j,k)  
               trk_trK = uppermet_uxx*kxx_p(i,j,k)+uppermet_uyy*kyy_p(i,j,k) 
     &              +uppermet_uzz*kzz_p(i,j,k)+2D0*(uppermet_uxy*kxy_p(i,j,k) 
     &              +uppermet_uxz*kxz_p(i,j,k)+uppermet_uyz*kyz_p(i,j,k))

               da_dxda = oo2dx*(alp(i+1,j,k)-alp(i-1,j,k))
               da_dyda = oo2dy*(alp(i,j+1,k)-alp(i,j-1,k))
               da_dzda = oo2dz*(alp(i,j,k+1)-alp(i,j,k-1))

               dda_dxxda = oodx2*(alp(i+1,j,k)-2D0*alp(i,j,k)+alp(i-1,j,k))
               dda_dyyda = oody2*(alp(i,j+1,k)-2D0*alp(i,j,k)+alp(i,j-1,k))
               dda_dzzda = oodz2*(alp(i,j,k+1)-2D0*alp(i,j,k)+alp(i,j,k-1))
               dda_dxyda = oo4dxdy*(alp(i+1,j+1,k)-alp(i+1,j-1,k)-alp(i-1,j+1,k)+alp(i-1,j-1,k))
               dda_dxzda = oo4dxdz*(alp(i+1,j,k+1)-alp(i+1,j,k-1)-alp(i-1,j,k+1)+alp(i-1,j,k-1))
               dda_dyzda = oo4dydz*(alp(i,j+1,k+1)-alp(i,j+1,k-1)-alp(i,j-1,k+1)+alp(i,j-1,k-1))

               cdcda_cdxxda = (dda_dxxda-gamma111*da_dxda-gamma211*da_dyda-gamma311*da_dzda)
               cdcda_cdxyda = (dda_dxyda-gamma112*da_dxda-gamma212*da_dyda-gamma312*da_dzda)
               cdcda_cdxzda = (dda_dxzda-gamma113*da_dxda-gamma213*da_dyda-gamma313*da_dzda)
               cdcda_cdyyda = (dda_dyyda-gamma122*da_dxda-gamma222*da_dyda-gamma322*da_dzda)
               cdcda_cdyzda = (dda_dyzda-gamma123*da_dxda-gamma223*da_dyda-gamma323*da_dzda)
               cdcda_cdzzda = (dda_dzzda-gamma133*da_dxda-gamma233*da_dyda-gamma333*da_dzda)

               adms_kxx(i,j,k) = alp(i,j,k)*(ricci_R11-2*KK11+kxx_p(i,j,k)*trk_trK)-cdcda_cdxxda 
               adms_kxy(i,j,k) = alp(i,j,k)*(ricci_R12-2*KK12+kxy_p(i,j,k)*trk_trK)-cdcda_cdxyda 
               adms_kxz(i,j,k) = alp(i,j,k)*(ricci_R13-2*KK13+kxz_p(i,j,k)*trk_trK)-cdcda_cdxzda 
               adms_kyy(i,j,k) = alp(i,j,k)*(ricci_R22-2*KK22+kyy_p(i,j,k)*trk_trK)-cdcda_cdyyda 
               adms_kyz(i,j,k) = alp(i,j,k)*(ricci_R23-2*KK23+kyz_p(i,j,k)*trk_trK)-cdcda_cdyzda 
               adms_kzz(i,j,k) = alp(i,j,k)*(ricci_R33-2*KK33+kzz_p(i,j,k)*trk_trK)-cdcda_cdzzda 

            end do
         end do
      end do

      return
      end