c DO NOT EDIT THIS FILE DIRECTLY! Rather edit adm_helper.pl which creates it
c Finite Difference the g and alp
c First derivatives
        dxalp = (alp(i+1,j,k) - alp(i-1,j,k))/(2.0D0*dx)
        dexgxx = (gxx(i+1,j,k) - gxx(i-1,j,k))/(2.0D0*dx)
        dexgxy = (gxy(i+1,j,k) - gxy(i-1,j,k))/(2.0D0*dx)
        dexgxz = (gxz(i+1,j,k) - gxz(i-1,j,k))/(2.0D0*dx)
        dexgyy = (gyy(i+1,j,k) - gyy(i-1,j,k))/(2.0D0*dx)
        dexgyz = (gyz(i+1,j,k) - gyz(i-1,j,k))/(2.0D0*dx)
        dexgzz = (gzz(i+1,j,k) - gzz(i-1,j,k))/(2.0D0*dx)
        dyalp = (alp(i,j+1,k) - alp(i,j-1,k))/(2.0D0*dy)
        deygxx = (gxx(i,j+1,k) - gxx(i,j-1,k))/(2.0D0*dy)
        deygxy = (gxy(i,j+1,k) - gxy(i,j-1,k))/(2.0D0*dy)
        deygxz = (gxz(i,j+1,k) - gxz(i,j-1,k))/(2.0D0*dy)
        deygyy = (gyy(i,j+1,k) - gyy(i,j-1,k))/(2.0D0*dy)
        deygyz = (gyz(i,j+1,k) - gyz(i,j-1,k))/(2.0D0*dy)
        deygzz = (gzz(i,j+1,k) - gzz(i,j-1,k))/(2.0D0*dy)
        dzalp = (alp(i,j,k+1) - alp(i,j,k-1))/(2.0D0*dz)
        dezgxx = (gxx(i,j,k+1) - gxx(i,j,k-1))/(2.0D0*dz)
        dezgxy = (gxy(i,j,k+1) - gxy(i,j,k-1))/(2.0D0*dz)
        dezgxz = (gxz(i,j,k+1) - gxz(i,j,k-1))/(2.0D0*dz)
        dezgyy = (gyy(i,j,k+1) - gyy(i,j,k-1))/(2.0D0*dz)
        dezgyz = (gyz(i,j,k+1) - gyz(i,j,k-1))/(2.0D0*dz)
        dezgzz = (gzz(i,j,k+1) - gzz(i,j,k-1))/(2.0D0*dz)
c Non-mixed second derivatives
        d2xxalp = (alp(i+1,j,k) + alp(i-1,j,k) - 2.0D0*alp(i,j,k))/(dx**2)
        d2xxgxx = (gxx(i+1,j,k) + gxx(i-1,j,k) - 2.0D0*gxx(i,j,k))/(dx**2)
        d2xxgxy = (gxy(i+1,j,k) + gxy(i-1,j,k) - 2.0D0*gxy(i,j,k))/(dx**2)
        d2xxgxz = (gxz(i+1,j,k) + gxz(i-1,j,k) - 2.0D0*gxz(i,j,k))/(dx**2)
        d2xxgyy = (gyy(i+1,j,k) + gyy(i-1,j,k) - 2.0D0*gyy(i,j,k))/(dx**2)
        d2xxgyz = (gyz(i+1,j,k) + gyz(i-1,j,k) - 2.0D0*gyz(i,j,k))/(dx**2)
        d2xxgzz = (gzz(i+1,j,k) + gzz(i-1,j,k) - 2.0D0*gzz(i,j,k))/(dx**2)
        d2yyalp = (alp(i,j+1,k) + alp(i,j-1,k) - 2.0D0*alp(i,j,k))/(dy**2)
        d2yygxx = (gxx(i,j+1,k) + gxx(i,j-1,k) - 2.0D0*gxx(i,j,k))/(dy**2)
        d2yygxy = (gxy(i,j+1,k) + gxy(i,j-1,k) - 2.0D0*gxy(i,j,k))/(dy**2)
        d2yygxz = (gxz(i,j+1,k) + gxz(i,j-1,k) - 2.0D0*gxz(i,j,k))/(dy**2)
        d2yygyy = (gyy(i,j+1,k) + gyy(i,j-1,k) - 2.0D0*gyy(i,j,k))/(dy**2)
        d2yygyz = (gyz(i,j+1,k) + gyz(i,j-1,k) - 2.0D0*gyz(i,j,k))/(dy**2)
        d2yygzz = (gzz(i,j+1,k) + gzz(i,j-1,k) - 2.0D0*gzz(i,j,k))/(dy**2)
        d2zzalp = (alp(i,j,k+1) + alp(i,j,k-1) - 2.0D0*alp(i,j,k))/(dz**2)
        d2zzgxx = (gxx(i,j,k+1) + gxx(i,j,k-1) - 2.0D0*gxx(i,j,k))/(dz**2)
        d2zzgxy = (gxy(i,j,k+1) + gxy(i,j,k-1) - 2.0D0*gxy(i,j,k))/(dz**2)
        d2zzgxz = (gxz(i,j,k+1) + gxz(i,j,k-1) - 2.0D0*gxz(i,j,k))/(dz**2)
        d2zzgyy = (gyy(i,j,k+1) + gyy(i,j,k-1) - 2.0D0*gyy(i,j,k))/(dz**2)
        d2zzgyz = (gyz(i,j,k+1) + gyz(i,j,k-1) - 2.0D0*gyz(i,j,k))/(dz**2)
        d2zzgzz = (gzz(i,j,k+1) + gzz(i,j,k-1) - 2.0D0*gzz(i,j,k))/(dz**2)
c Mixed second derivatives 
        d2xyalp = (alp(i+1,j+1,k) + alp(i-1,j-1,k) - 
     $  alp(i-1,j+1,k) - alp(i+1,j-1,k))/((4.0D0*dx*dy))
        d2xygxx = (gxx(i+1,j+1,k) + gxx(i-1,j-1,k) - 
     $  gxx(i-1,j+1,k) - gxx(i+1,j-1,k))/((4.0D0*dx*dy))
        d2xygxy = (gxy(i+1,j+1,k) + gxy(i-1,j-1,k) - 
     $  gxy(i-1,j+1,k) - gxy(i+1,j-1,k))/((4.0D0*dx*dy))
        d2xygxz = (gxz(i+1,j+1,k) + gxz(i-1,j-1,k) - 
     $  gxz(i-1,j+1,k) - gxz(i+1,j-1,k))/((4.0D0*dx*dy))
        d2xygyy = (gyy(i+1,j+1,k) + gyy(i-1,j-1,k) - 
     $  gyy(i-1,j+1,k) - gyy(i+1,j-1,k))/((4.0D0*dx*dy))
        d2xygyz = (gyz(i+1,j+1,k) + gyz(i-1,j-1,k) - 
     $  gyz(i-1,j+1,k) - gyz(i+1,j-1,k))/((4.0D0*dx*dy))
        d2xygzz = (gzz(i+1,j+1,k) + gzz(i-1,j-1,k) - 
     $  gzz(i-1,j+1,k) - gzz(i+1,j-1,k))/((4.0D0*dx*dy))
        d2xzalp = (alp(i+1,j,k+1) + alp(i-1,j,k-1) - 
     $  alp(i-1,j,k+1) - alp(i+1,j,k-1))/((4.0D0*dx*dz))
        d2xzgxx = (gxx(i+1,j,k+1) + gxx(i-1,j,k-1) - 
     $  gxx(i-1,j,k+1) - gxx(i+1,j,k-1))/((4.0D0*dx*dz))
        d2xzgxy = (gxy(i+1,j,k+1) + gxy(i-1,j,k-1) - 
     $  gxy(i-1,j,k+1) - gxy(i+1,j,k-1))/((4.0D0*dx*dz))
        d2xzgxz = (gxz(i+1,j,k+1) + gxz(i-1,j,k-1) - 
     $  gxz(i-1,j,k+1) - gxz(i+1,j,k-1))/((4.0D0*dx*dz))
        d2xzgyy = (gyy(i+1,j,k+1) + gyy(i-1,j,k-1) - 
     $  gyy(i-1,j,k+1) - gyy(i+1,j,k-1))/((4.0D0*dx*dz))
        d2xzgyz = (gyz(i+1,j,k+1) + gyz(i-1,j,k-1) - 
     $  gyz(i-1,j,k+1) - gyz(i+1,j,k-1))/((4.0D0*dx*dz))
        d2xzgzz = (gzz(i+1,j,k+1) + gzz(i-1,j,k-1) - 
     $  gzz(i-1,j,k+1) - gzz(i+1,j,k-1))/((4.0D0*dx*dz))
        d2yzalp = (alp(i,j+1,k+1) + alp(i,j-1,k-1) - 
     $  alp(i,j-1,k+1) - alp(i,j+1,k-1))/((4.0D0*dy*dz))
        d2yzgxx = (gxx(i,j+1,k+1) + gxx(i,j-1,k-1) - 
     $  gxx(i,j-1,k+1) - gxx(i,j+1,k-1))/((4.0D0*dy*dz))
        d2yzgxy = (gxy(i,j+1,k+1) + gxy(i,j-1,k-1) - 
     $  gxy(i,j-1,k+1) - gxy(i,j+1,k-1))/((4.0D0*dy*dz))
        d2yzgxz = (gxz(i,j+1,k+1) + gxz(i,j-1,k-1) - 
     $  gxz(i,j-1,k+1) - gxz(i,j+1,k-1))/((4.0D0*dy*dz))
        d2yzgyy = (gyy(i,j+1,k+1) + gyy(i,j-1,k-1) - 
     $  gyy(i,j-1,k+1) - gyy(i,j+1,k-1))/((4.0D0*dy*dz))
        d2yzgyz = (gyz(i,j+1,k+1) + gyz(i,j-1,k-1) - 
     $  gyz(i,j-1,k+1) - gyz(i,j+1,k-1))/((4.0D0*dy*dz))
        d2yzgzz = (gzz(i,j+1,k+1) + gzz(i,j-1,k-1) - 
     $  gzz(i,j-1,k+1) - gzz(i,j+1,k-1))/((4.0D0*dy*dz))