


























































































































 


































































































































































 
































































































































































































 



















































 










 








































































































































































































































 










                                
                                

                                
                                

                                
                                

                                



















































































subroutine NP_calcPsiGF( cctk_dim,cctk_gsh,cctk_lsh,cctk_lbnd,cctk_ubnd,cctk_tile_min,cctk_tile_max,cctk_ash,cctk_alignment,cctk_al&
  &ignment_offset,cctk_from,cctk_to,cctk_bbox,cctk_delta_time,cctk_time,cctk_delta_space,cctk_origin_space,cctk_levfac,cctk_levoff,&
  &cctk_levoffdenom,cctk_timefac,cctk_convlevel,cctk_convfac,cctk_nghostzones,cctk_iteration,cctkGH,cctk_ash1,cctk_ash2,cctk_ash3, &
  &psi4im,psi4im_p,psi4im_p_p,psi4re,psi4re_p,psi4re_p_p, alp,alp_p,alp_p_p,betax,betax_p,betax_p_p,betay,betay_p,betay_p_p,betaz,b&
  &etaz_p,betaz_p_p,coarse_dx,coarse_dy,coarse_dz,dtalp,dtalp_p,dtalp_p_p,dtbetax,dtbetax_p,dtbetax_p_p,dtbetay,dtbetay_p,dtbetay_p&
  &_p,dtbetaz,dtbetaz_p,dtbetaz_p_p,dtlapse_state,dtshift_state,gxx,gxx_p,gxx_p_p,gxy,gxy_p,gxy_p_p,gxz,gxz_p,gxz_p_p,gyy,gyy_p,gyy&
  &_p_p,gyz,gyz_p,gyz_p_p,gzz,gzz_p,gzz_p_p,kxx,kxx_p,kxx_p_p,kxy,kxy_p,kxy_p_p,kxz,kxz_p,kxz_p_p,kyy,kyy_p,kyy_p_p,kyz,kyz_p,kyz_p&
  &_p,kzz,kzz_p,kzz_p_p,r,shift_state,x,y,z )

  implicit none

  INTEGER cctk_dim
  integer, parameter :: cctki_use_cctk_dim = kind(cctk_dim)
  INTEGER cctk_gsh (cctk_dim)
  integer, parameter :: cctki_use_cctk_gsh = kind(cctk_gsh)
  INTEGER cctk_lsh (cctk_dim)
  integer, parameter :: cctki_use_cctk_lsh = kind(cctk_lsh)
  INTEGER cctk_lbnd (cctk_dim)
  integer, parameter :: cctki_use_cctk_lbnd = kind(cctk_lbnd)
  INTEGER cctk_ubnd (cctk_dim)
  integer, parameter :: cctki_use_cctk_ubnd = kind(cctk_ubnd)
  INTEGER cctk_tile_min (cctk_dim)
  integer, parameter :: cctki_use_cctk_tile_min = kind(cctk_tile_min)
  INTEGER cctk_tile_max (cctk_dim)
  integer, parameter :: cctki_use_cctk_tile_max = kind(cctk_tile_max)
  INTEGER cctk_ash (cctk_dim)
  integer, parameter :: cctki_use_cctk_ash = kind(cctk_ash)
  INTEGER cctk_alignment
  integer, parameter :: cctki_use_cctk_alignment = kind(cctk_alignment)
  INTEGER cctk_alignment_offset
  integer, parameter :: cctki_use_cctk_alignment_offset = kind(cctk_alignment_offset)
  INTEGER cctk_from (cctk_dim)
  integer, parameter :: cctki_use_cctk_from = kind(cctk_from)
  INTEGER cctk_to (cctk_dim)
  integer, parameter :: cctki_use_cctk_to = kind(cctk_to)
  INTEGER cctk_bbox (2*cctk_dim)
  integer, parameter :: cctki_use_cctk_bbox = kind(cctk_bbox)
  REAL*8 cctk_delta_time
  integer, parameter :: cctki_use_cctk_delta_time = kind(cctk_delta_time)
  REAL*8 cctk_time
  integer, parameter :: cctki_use_cctk_time = kind(cctk_time)
  REAL*8 cctk_delta_space (cctk_dim)
  integer, parameter :: cctki_use_cctk_delta_space = kind(cctk_delta_space)
  REAL*8 cctk_origin_space (cctk_dim)
  integer, parameter :: cctki_use_cctk_origin_space = kind(cctk_origin_space)
  INTEGER cctk_levfac (cctk_dim)
  integer, parameter :: cctki_use_cctk_levfac = kind(cctk_levfac)
  INTEGER cctk_levoff (cctk_dim)
  integer, parameter :: cctki_use_cctk_levoff = kind(cctk_levoff)
  INTEGER cctk_levoffdenom (cctk_dim)
  integer, parameter :: cctki_use_cctk_levoffdenom = kind(cctk_levoffdenom)
  INTEGER cctk_timefac
  integer, parameter :: cctki_use_cctk_timefac = kind(cctk_timefac)
  INTEGER cctk_convlevel
  integer, parameter :: cctki_use_cctk_convlevel = kind(cctk_convlevel)
  INTEGER cctk_convfac
  integer, parameter :: cctki_use_cctk_convfac = kind(cctk_convfac)
  INTEGER cctk_nghostzones (cctk_dim)
  integer, parameter :: cctki_use_cctk_nghostzones = kind(cctk_nghostzones)
  INTEGER cctk_iteration
  integer, parameter :: cctki_use_cctk_iteration = kind(cctk_iteration)
  integer*8 cctkGH
  integer, parameter :: cctki_use_cctkGH = kind(cctkGH)
  INTEGER cctk_ash1
  integer, parameter :: cctki_use_cctk_ash1 = kind(cctk_ash1)
  INTEGER cctk_ash2
  integer, parameter :: cctki_use_cctk_ash2 = kind(cctk_ash2)
  INTEGER cctk_ash3
  integer, parameter :: cctki_use_cctk_ash3 = kind(cctk_ash3)
  REAL*8 psi4im (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_psi4im = kind(psi4im)
  REAL*8 psi4im_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_psi4im_p = kind(psi4im_p)
  REAL*8 psi4im_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_psi4im_p_p = kind(psi4im_p_p)
  REAL*8 psi4re (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_psi4re = kind(psi4re)
  REAL*8 psi4re_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_psi4re_p = kind(psi4re_p)
  REAL*8 psi4re_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_psi4re_p_p = kind(psi4re_p_p)
  REAL*8 alp (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_alp = kind(alp)
  REAL*8 alp_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_alp_p = kind(alp_p)
  REAL*8 alp_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_alp_p_p = kind(alp_p_p)
  REAL*8 betax (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_betax = kind(betax)
  REAL*8 betax_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_betax_p = kind(betax_p)
  REAL*8 betax_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_betax_p_p = kind(betax_p_p)
  REAL*8 betay (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_betay = kind(betay)
  REAL*8 betay_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_betay_p = kind(betay_p)
  REAL*8 betay_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_betay_p_p = kind(betay_p_p)
  REAL*8 betaz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_betaz = kind(betaz)
  REAL*8 betaz_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_betaz_p = kind(betaz_p)
  REAL*8 betaz_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_betaz_p_p = kind(betaz_p_p)
  REAL*8 coarse_dx
  integer, parameter :: cctki_use_coarse_dx = kind(coarse_dx)
  REAL*8 coarse_dy
  integer, parameter :: cctki_use_coarse_dy = kind(coarse_dy)
  REAL*8 coarse_dz
  integer, parameter :: cctki_use_coarse_dz = kind(coarse_dz)
  REAL*8 dtalp (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_dtalp = kind(dtalp)
  REAL*8 dtalp_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_dtalp_p = kind(dtalp_p)
  REAL*8 dtalp_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_dtalp_p_p = kind(dtalp_p_p)
  REAL*8 dtbetax (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_dtbetax = kind(dtbetax)
  REAL*8 dtbetax_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_dtbetax_p = kind(dtbetax_p)
  REAL*8 dtbetax_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_dtbetax_p_p = kind(dtbetax_p_p)
  REAL*8 dtbetay (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_dtbetay = kind(dtbetay)
  REAL*8 dtbetay_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_dtbetay_p = kind(dtbetay_p)
  REAL*8 dtbetay_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_dtbetay_p_p = kind(dtbetay_p_p)
  REAL*8 dtbetaz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_dtbetaz = kind(dtbetaz)
  REAL*8 dtbetaz_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_dtbetaz_p = kind(dtbetaz_p)
  REAL*8 dtbetaz_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_dtbetaz_p_p = kind(dtbetaz_p_p)
  INTEGER*4 dtlapse_state
  integer, parameter :: cctki_use_dtlapse_state = kind(dtlapse_state)
  INTEGER*4 dtshift_state
  integer, parameter :: cctki_use_dtshift_state = kind(dtshift_state)
  REAL*8 gxx (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_gxx = kind(gxx)
  REAL*8 gxx_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_gxx_p = kind(gxx_p)
  REAL*8 gxx_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_gxx_p_p = kind(gxx_p_p)
  REAL*8 gxy (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_gxy = kind(gxy)
  REAL*8 gxy_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_gxy_p = kind(gxy_p)
  REAL*8 gxy_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_gxy_p_p = kind(gxy_p_p)
  REAL*8 gxz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_gxz = kind(gxz)
  REAL*8 gxz_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_gxz_p = kind(gxz_p)
  REAL*8 gxz_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_gxz_p_p = kind(gxz_p_p)
  REAL*8 gyy (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_gyy = kind(gyy)
  REAL*8 gyy_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_gyy_p = kind(gyy_p)
  REAL*8 gyy_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_gyy_p_p = kind(gyy_p_p)
  REAL*8 gyz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_gyz = kind(gyz)
  REAL*8 gyz_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_gyz_p = kind(gyz_p)
  REAL*8 gyz_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_gyz_p_p = kind(gyz_p_p)
  REAL*8 gzz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_gzz = kind(gzz)
  REAL*8 gzz_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_gzz_p = kind(gzz_p)
  REAL*8 gzz_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_gzz_p_p = kind(gzz_p_p)
  REAL*8 kxx (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_kxx = kind(kxx)
  REAL*8 kxx_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_kxx_p = kind(kxx_p)
  REAL*8 kxx_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_kxx_p_p = kind(kxx_p_p)
  REAL*8 kxy (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_kxy = kind(kxy)
  REAL*8 kxy_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_kxy_p = kind(kxy_p)
  REAL*8 kxy_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_kxy_p_p = kind(kxy_p_p)
  REAL*8 kxz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_kxz = kind(kxz)
  REAL*8 kxz_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_kxz_p = kind(kxz_p)
  REAL*8 kxz_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_kxz_p_p = kind(kxz_p_p)
  REAL*8 kyy (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_kyy = kind(kyy)
  REAL*8 kyy_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_kyy_p = kind(kyy_p)
  REAL*8 kyy_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_kyy_p_p = kind(kyy_p_p)
  REAL*8 kyz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_kyz = kind(kyz)
  REAL*8 kyz_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_kyz_p = kind(kyz_p)
  REAL*8 kyz_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_kyz_p_p = kind(kyz_p_p)
  REAL*8 kzz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_kzz = kind(kzz)
  REAL*8 kzz_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_kzz_p = kind(kzz_p)
  REAL*8 kzz_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_kzz_p_p = kind(kzz_p_p)
  REAL*8 r (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_r = kind(r)
  INTEGER*4 shift_state
  integer, parameter :: cctki_use_shift_state = kind(shift_state)
  REAL*8 x (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_x = kind(x)
  REAL*8 y (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_y = kind(y)
  REAL*8 z (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_z = kind(z)
  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
  interface
  INTEGER*4 function Boundary_SelectGroupForBC (GH, faces, boundary_width, table_handle, group_name, bc_name)
  implicit none
  integer*8 GH
  INTEGER*4 faces
  INTEGER*4 boundary_width
  INTEGER*4 table_handle
  character(*) group_name
  character(*) bc_name
  end function Boundary_SelectGroupForBC
  end interface
  INTEGER*4 NP_order
  integer, parameter :: cctki_use_NP_order = kind(NP_order)
  INTEGER*4 calculate_NP_every
  integer, parameter :: cctki_use_calculate_NP_every = kind(calculate_NP_every)
  COMMON /NPScalarspriv/NP_order, calculate_NP_every

  REAL*8 gd(3,3), gu(3,3), kd(3,3), detgd
  REAL*8 d1_gd(3,3,3), d1_gu(3,3,3), d2_gd(3,3,3,3), d1_kd(3,3,3),      &
            cd1_kd(3,3,3)
  REAL*8 cf1(3,3,3), cf2(3,3,3), ri(3,3), d1_cf2(3,3,3,3)
  REAL*8 dx(3), xx(3)
  REAL*8 dx12, dy12, dz12, dxsq12, dysq12, dzsq12, dxdy144, dxdz144,    &
            dydz144
  REAL*8 odx60, ody60, odz60, odxsq180, odysq180, odzsq180,&
            odxdy3600, odxdz3600, odydz3600
  REAL*8 dx2, dy2, dz2, dxsq, dysq, dzsq, dxdy4, dxdz4, dydz4
  REAL*8 u_vec(3), v_vec(3), w_vec(3), ud_vec(3), dotp1, dotp2
  REAL*8 eps_lc_d(3,3,3), eps_lc_u(3,3,3), elec(3,3), mag(3,3), electr, magtr

  INTEGER*4  i, j, k, m, n, p, a, b, c, d

  if (calculate_NP_every .le. 0) then
     return
  end if

  if (MOD(cctk_iteration, calculate_NP_every) .ne. 0 ) then
     return
  endif

  dx(:)   = (cctk_delta_space(:)/cctk_levfac(:))

  
  odx60 = 1 / (60 * (cctk_delta_space(1)/cctk_levfac(1)))
  ody60 = 1 / (60 * (cctk_delta_space(2)/cctk_levfac(2)))
  odz60 = 1 / (60 * (cctk_delta_space(3)/cctk_levfac(3)))

  odxsq180 = 1 / (180*(cctk_delta_space(1)/cctk_levfac(1))**2)
  odysq180 = 1 / (180*(cctk_delta_space(2)/cctk_levfac(2))**2)
  odzsq180 = 1 / (180*(cctk_delta_space(3)/cctk_levfac(3))**2)

  odxdy3600 = 1 / (3600*(cctk_delta_space(1)/cctk_levfac(1))*(cctk_delta_space(2)/cctk_levfac(2)))
  odxdz3600 = 1 / (3600*(cctk_delta_space(1)/cctk_levfac(1))*(cctk_delta_space(3)/cctk_levfac(3)))
  odydz3600 = 1 / (3600*(cctk_delta_space(2)/cctk_levfac(2))*(cctk_delta_space(3)/cctk_levfac(3)))

  
  dx12    = 12 * dx(1)
  dy12    = 12 * dx(2)
  dz12    = 12 * dx(3)
  dxsq12  = 12 * dx(1)**2
  dysq12  = 12 * dx(2)**2
  dzsq12  = 12 * dx(3)**2
  dxdy144 = dx12 * dy12
  dxdz144 = dx12 * dz12
  dydz144 = dy12 * dz12

  
  dx2   = 2 * dx(1)
  dy2   = 2 * dx(2)
  dz2   = 2 * dx(3)
  dxsq  = dx(1)**2
  dysq  = dx(2)**2
  dzsq  = dx(3)**2
  dxdy4 = 4*dx(1) * dx(2)
  dxdz4 = 4*dx(1) * dx(3)
  dydz4 = 4*dx(2) * dx(3)

  
  psi4re = 0
  psi4im = 0

  !$OMP PARALLEL DO COLLAPSE(3) &
  !$OMP PRIVATE(gd, kd, &
  !$OMP detgd, gu, &
  !$OMP d1_gd, d1_kd, d2_gd, &
  !$OMP cf1, cf2, cd1_kd, d1_gu, d1_cf2, &
  !$OMP ri, eps_lc_u, eps_lc_d, &
  !$OMP xx, u_vec, v_vec, w_vec, dotp1, dotp2, &
  !$OMP ud_vec, elec, mag, electr, magtr, &
  !$OMP i, j, k, &
  !$OMP a, b, c, m, n)
  do k = 1+cctk_nghostzones(3), cctk_lsh(3)-cctk_nghostzones(3)
  do j = 1+cctk_nghostzones(2), cctk_lsh(2)-cctk_nghostzones(2)
  do i = 1+cctk_nghostzones(1), cctk_lsh(1)-cctk_nghostzones(1)

    
    gd(1,1) = gxx(i,j,k)
    gd(1,2) = gxy(i,j,k)
    gd(1,3) = gxz(i,j,k)
    gd(2,2) = gyy(i,j,k)
    gd(2,3) = gyz(i,j,k)
    gd(3,3) = gzz(i,j,k)
    gd(2,1) = gd(1,2)
    gd(3,1) = gd(1,3)
    gd(3,2) = gd(2,3)

    kd(1,1) = kxx(i,j,k)
    kd(1,2) = kxy(i,j,k)
    kd(1,3) = kxz(i,j,k)
    kd(2,2) = kyy(i,j,k)
    kd(2,3) = kyz(i,j,k)
    kd(3,3) = kzz(i,j,k)
    kd(2,1) = kd(1,2)
    kd(3,1) = kd(1,3)
    kd(3,2) = kd(2,3)
    


    
    detgd =       gd(1,1) * gd(2,2) * gd(3,3)                                &
            + 2 * gd(1,2) * gd(1,3) * gd(2,3)                                &
            -     gd(1,1) * gd(2,3) ** 2                                     &
            -     gd(2,2) * gd(1,3) ** 2                                     &
            -     gd(3,3) * gd(1,2) ** 2
    gu(1,1) = (gd(2,2) * gd(3,3) - gd(2,3) ** 2     ) / detgd
    gu(2,2) = (gd(1,1) * gd(3,3) - gd(1,3) ** 2     ) / detgd
    gu(3,3) = (gd(1,1) * gd(2,2) - gd(1,2) ** 2     ) / detgd
    gu(1,2) = (gd(1,3) * gd(2,3) - gd(1,2) * gd(3,3)) / detgd
    gu(1,3) = (gd(1,2) * gd(2,3) - gd(1,3) * gd(2,2)) / detgd
    gu(2,3) = (gd(1,3) * gd(1,2) - gd(2,3) * gd(1,1)) / detgd
    gu(2,1) = gu(1,2)
    gu(3,1) = gu(1,3)
    gu(3,2) = gu(2,3)
    


    
    if ( NP_order == 6 ) then
    
    
    
      d1_gd(1,1,1) = (  gxx(i+3,j,k) - 9*gxx(i+2,j,k) + 45*gxx(i+1,j,k)      &
                      - gxx(i-3,j,k) + 9*gxx(i-2,j,k) - 45*gxx(i-1,j,k) ) * odx60
      d1_gd(1,2,1) = (  gxy(i+3,j,k) - 9*gxy(i+2,j,k) + 45*gxy(i+1,j,k)      &
                      - gxy(i-3,j,k) + 9*gxy(i-2,j,k) - 45*gxy(i-1,j,k) ) * odx60
      d1_gd(1,3,1) = (  gxz(i+3,j,k) - 9*gxz(i+2,j,k) + 45*gxz(i+1,j,k)      &
                      - gxz(i-3,j,k) + 9*gxz(i-2,j,k) - 45*gxz(i-1,j,k) ) * odx60
      d1_gd(2,2,1) = (  gyy(i+3,j,k) - 9*gyy(i+2,j,k) + 45*gyy(i+1,j,k)      &
                      - gyy(i-3,j,k) + 9*gyy(i-2,j,k) - 45*gyy(i-1,j,k) ) * odx60
      d1_gd(2,3,1) = (  gyz(i+3,j,k) - 9*gyz(i+2,j,k) + 45*gyz(i+1,j,k)      &
                      - gyz(i-3,j,k) + 9*gyz(i-2,j,k) - 45*gyz(i-1,j,k) ) * odx60
      d1_gd(3,3,1) = (  gzz(i+3,j,k) - 9*gzz(i+2,j,k) + 45*gzz(i+1,j,k)      &
                      - gzz(i-3,j,k) + 9*gzz(i-2,j,k) - 45*gzz(i-1,j,k) ) * odx60

      d1_gd(1,1,2) = (  gxx(i,j+3,k) - 9*gxx(i,j+2,k) + 45*gxx(i,j+1,k)      &
                      - gxx(i,j-3,k) + 9*gxx(i,j-2,k) - 45*gxx(i,j-1,k) ) * ody60
      d1_gd(1,2,2) = (  gxy(i,j+3,k) - 9*gxy(i,j+2,k) + 45*gxy(i,j+1,k)      &
                      - gxy(i,j-3,k) + 9*gxy(i,j-2,k) - 45*gxy(i,j-1,k) ) * ody60
      d1_gd(1,3,2) = (  gxz(i,j+3,k) - 9*gxz(i,j+2,k) + 45*gxz(i,j+1,k)      &
                      - gxz(i,j-3,k) + 9*gxz(i,j-2,k) - 45*gxz(i,j-1,k) ) * ody60
      d1_gd(2,2,2) = (  gyy(i,j+3,k) - 9*gyy(i,j+2,k) + 45*gyy(i,j+1,k)      &
                      - gyy(i,j-3,k) + 9*gyy(i,j-2,k) - 45*gyy(i,j-1,k) ) * ody60
      d1_gd(2,3,2) = (  gyz(i,j+3,k) - 9*gyz(i,j+2,k) + 45*gyz(i,j+1,k)      &
                      - gyz(i,j-3,k) + 9*gyz(i,j-2,k) - 45*gyz(i,j-1,k) ) * ody60
      d1_gd(3,3,2) = (  gzz(i,j+3,k) - 9*gzz(i,j+2,k) + 45*gzz(i,j+1,k)      &
                      - gzz(i,j-3,k) + 9*gzz(i,j-2,k) - 45*gzz(i,j-1,k) ) * ody60

      d1_gd(1,1,3) = (  gxx(i,j,k+3) - 9*gxx(i,j,k+2) + 45*gxx(i,j,k+1)      &
                      - gxx(i,j,k-3) + 9*gxx(i,j,k-2) - 45*gxx(i,j,k-1) ) * odz60
      d1_gd(1,2,3) = (  gxy(i,j,k+3) - 9*gxy(i,j,k+2) + 45*gxy(i,j,k+1)      &
                      - gxy(i,j,k-3) + 9*gxy(i,j,k-2) - 45*gxy(i,j,k-1) ) * odz60
      d1_gd(1,3,3) = (  gxz(i,j,k+3) - 9*gxz(i,j,k+2) + 45*gxz(i,j,k+1)      &
                      - gxz(i,j,k-3) + 9*gxz(i,j,k-2) - 45*gxz(i,j,k-1) ) * odz60
      d1_gd(2,2,3) = (  gyy(i,j,k+3) - 9*gyy(i,j,k+2) + 45*gyy(i,j,k+1)      &
                      - gyy(i,j,k-3) + 9*gyy(i,j,k-2) - 45*gyy(i,j,k-1) ) * odz60
      d1_gd(2,3,3) = (  gyz(i,j,k+3) - 9*gyz(i,j,k+2) + 45*gyz(i,j,k+1)      &
                      - gyz(i,j,k-3) + 9*gyz(i,j,k-2) - 45*gyz(i,j,k-1) ) * odz60
      d1_gd(3,3,3) = (  gzz(i,j,k+3) - 9*gzz(i,j,k+2) + 45*gzz(i,j,k+1)      &
                      - gzz(i,j,k-3) + 9*gzz(i,j,k-2) - 45*gzz(i,j,k-1) ) * odz60

      d1_gd(2,1,:) = d1_gd(1,2,:)
      d1_gd(3,1,:) = d1_gd(1,3,:)
      d1_gd(3,2,:) = d1_gd(2,3,:)


    
      d1_kd(1,1,1) = (  kxx(i+3,j,k) - 9*kxx(i+2,j,k) + 45*kxx(i+1,j,k)      &
                      - kxx(i-3,j,k) + 9*kxx(i-2,j,k) - 45*kxx(i-1,j,k) ) * odx60
      d1_kd(1,2,1) = (  kxy(i+3,j,k) - 9*kxy(i+2,j,k) + 45*kxy(i+1,j,k)      &
                      - kxy(i-3,j,k) + 9*kxy(i-2,j,k) - 45*kxy(i-1,j,k) ) * odx60
      d1_kd(1,3,1) = (  kxz(i+3,j,k) - 9*kxz(i+2,j,k) + 45*kxz(i+1,j,k)      &
                      - kxz(i-3,j,k) + 9*kxz(i-2,j,k) - 45*kxz(i-1,j,k) ) * odx60
      d1_kd(2,2,1) = (  kyy(i+3,j,k) - 9*kyy(i+2,j,k) + 45*kyy(i+1,j,k)      &
                      - kyy(i-3,j,k) + 9*kyy(i-2,j,k) - 45*kyy(i-1,j,k) ) * odx60
      d1_kd(2,3,1) = (  kyz(i+3,j,k) - 9*kyz(i+2,j,k) + 45*kyz(i+1,j,k)      &
                      - kyz(i-3,j,k) + 9*kyz(i-2,j,k) - 45*kyz(i-1,j,k) ) * odx60
      d1_kd(3,3,1) = (  kzz(i+3,j,k) - 9*kzz(i+2,j,k) + 45*kzz(i+1,j,k)      &
                      - kzz(i-3,j,k) + 9*kzz(i-2,j,k) - 45*kzz(i-1,j,k) ) * odx60

      d1_kd(1,1,2) = (  kxx(i,j+3,k) - 9*kxx(i,j+2,k) + 45*kxx(i,j+1,k)      &
                      - kxx(i,j-3,k) + 9*kxx(i,j-2,k) - 45*kxx(i,j-1,k) ) * ody60
      d1_kd(1,2,2) = (  kxy(i,j+3,k) - 9*kxy(i,j+2,k) + 45*kxy(i,j+1,k)      &
                      - kxy(i,j-3,k) + 9*kxy(i,j-2,k) - 45*kxy(i,j-1,k) ) * ody60
      d1_kd(1,3,2) = (  kxz(i,j+3,k) - 9*kxz(i,j+2,k) + 45*kxz(i,j+1,k)      &
                      - kxz(i,j-3,k) + 9*kxz(i,j-2,k) - 45*kxz(i,j-1,k) ) * ody60
      d1_kd(2,2,2) = (  kyy(i,j+3,k) - 9*kyy(i,j+2,k) + 45*kyy(i,j+1,k)      &
                      - kyy(i,j-3,k) + 9*kyy(i,j-2,k) - 45*kyy(i,j-1,k) ) * ody60
      d1_kd(2,3,2) = (  kyz(i,j+3,k) - 9*kyz(i,j+2,k) + 45*kyz(i,j+1,k)      &
                      - kyz(i,j-3,k) + 9*kyz(i,j-2,k) - 45*kyz(i,j-1,k) ) * ody60
      d1_kd(3,3,2) = (  kzz(i,j+3,k) - 9*kzz(i,j+2,k) + 45*kzz(i,j+1,k)      &
                      - kzz(i,j-3,k) + 9*kzz(i,j-2,k) - 45*kzz(i,j-1,k) ) * ody60

      d1_kd(1,1,3) = (  kxx(i,j,k+3) - 9*kxx(i,j,k+2) + 45*kxx(i,j,k+1)      &
                      - kxx(i,j,k-3) + 9*kxx(i,j,k-2) - 45*kxx(i,j,k-1) ) * odz60
      d1_kd(1,2,3) = (  kxy(i,j,k+3) - 9*kxy(i,j,k+2) + 45*kxy(i,j,k+1)      &
                      - kxy(i,j,k-3) + 9*kxy(i,j,k-2) - 45*kxy(i,j,k-1) ) * odz60
      d1_kd(1,3,3) = (  kxz(i,j,k+3) - 9*kxz(i,j,k+2) + 45*kxz(i,j,k+1)      &
                      - kxz(i,j,k-3) + 9*kxz(i,j,k-2) - 45*kxz(i,j,k-1) ) * odz60
      d1_kd(2,2,3) = (  kyy(i,j,k+3) - 9*kyy(i,j,k+2) + 45*kyy(i,j,k+1)      &
                      - kyy(i,j,k-3) + 9*kyy(i,j,k-2) - 45*kyy(i,j,k-1) ) * odz60
      d1_kd(2,3,3) = (  kyz(i,j,k+3) - 9*kyz(i,j,k+2) + 45*kyz(i,j,k+1)      &
                      - kyz(i,j,k-3) + 9*kyz(i,j,k-2) - 45*kyz(i,j,k-1) ) * odz60
      d1_kd(3,3,3) = (  kzz(i,j,k+3) - 9*kzz(i,j,k+2) + 45*kzz(i,j,k+1)      &
                      - kzz(i,j,k-3) + 9*kzz(i,j,k-2) - 45*kzz(i,j,k-1) ) * odz60

      d1_kd(2,1,:) = d1_kd(1,2,:)
      d1_kd(3,1,:) = d1_kd(1,3,:)
      d1_kd(3,2,:) = d1_kd(2,3,:)



    
    
      d2_gd(1,1,1,1) = (  2*gxx(i+3,j,k) - 27*gxx(i+2,j,k) + 270*gxx(i+1,j,k) - 490*gxx(i,j,k)&
                        + 2*gxx(i-3,j,k) - 27*gxx(i-2,j,k) + 270*gxx(i-1,j,k) ) * odxsq180
      d2_gd(1,2,1,1) = (  2*gxy(i+3,j,k) - 27*gxy(i+2,j,k) + 270*gxy(i+1,j,k) - 490*gxy(i,j,k)&
                        + 2*gxy(i-3,j,k) - 27*gxy(i-2,j,k) + 270*gxy(i-1,j,k) ) * odxsq180
      d2_gd(1,3,1,1) = (  2*gxz(i+3,j,k) - 27*gxz(i+2,j,k) + 270*gxz(i+1,j,k) - 490*gxz(i,j,k)&
                        + 2*gxz(i-3,j,k) - 27*gxz(i-2,j,k) + 270*gxz(i-1,j,k) ) * odxsq180
      d2_gd(2,2,1,1) = (  2*gyy(i+3,j,k) - 27*gyy(i+2,j,k) + 270*gyy(i+1,j,k) - 490*gyy(i,j,k)&
                        + 2*gyy(i-3,j,k) - 27*gyy(i-2,j,k) + 270*gyy(i-1,j,k) ) * odxsq180
      d2_gd(2,3,1,1) = (  2*gyz(i+3,j,k) - 27*gyz(i+2,j,k) + 270*gyz(i+1,j,k) - 490*gyz(i,j,k)&
                        + 2*gyz(i-3,j,k) - 27*gyz(i-2,j,k) + 270*gyz(i-1,j,k) ) * odxsq180
      d2_gd(3,3,1,1) = (  2*gzz(i+3,j,k) - 27*gzz(i+2,j,k) + 270*gzz(i+1,j,k) - 490*gzz(i,j,k)&
                        + 2*gzz(i-3,j,k) - 27*gzz(i-2,j,k) + 270*gzz(i-1,j,k) ) * odxsq180

      d2_gd(1,1,2,2) = (  2*gxx(i,j+3,k) - 27*gxx(i,j+2,k) + 270*gxx(i,j+1,k) - 490*gxx(i,j,k)&
                        + 2*gxx(i,j-3,k) - 27*gxx(i,j-2,k) + 270*gxx(i,j-1,k) ) * odysq180
      d2_gd(1,2,2,2) = (  2*gxy(i,j+3,k) - 27*gxy(i,j+2,k) + 270*gxy(i,j+1,k) - 490*gxy(i,j,k)&
                        + 2*gxy(i,j-3,k) - 27*gxy(i,j-2,k) + 270*gxy(i,j-1,k) ) * odysq180
      d2_gd(1,3,2,2) = (  2*gxz(i,j+3,k) - 27*gxz(i,j+2,k) + 270*gxz(i,j+1,k) - 490*gxz(i,j,k)&
                        + 2*gxz(i,j-3,k) - 27*gxz(i,j-2,k) + 270*gxz(i,j-1,k) ) * odysq180
      d2_gd(2,2,2,2) = (  2*gyy(i,j+3,k) - 27*gyy(i,j+2,k) + 270*gyy(i,j+1,k) - 490*gyy(i,j,k)&
                        + 2*gyy(i,j-3,k) - 27*gyy(i,j-2,k) + 270*gyy(i,j-1,k) ) * odysq180
      d2_gd(2,3,2,2) = (  2*gyz(i,j+3,k) - 27*gyz(i,j+2,k) + 270*gyz(i,j+1,k) - 490*gyz(i,j,k)&
                        + 2*gyz(i,j-3,k) - 27*gyz(i,j-2,k) + 270*gyz(i,j-1,k) ) * odysq180
      d2_gd(3,3,2,2) = (  2*gzz(i,j+3,k) - 27*gzz(i,j+2,k) + 270*gzz(i,j+1,k) - 490*gzz(i,j,k)&
                        + 2*gzz(i,j-3,k) - 27*gzz(i,j-2,k) + 270*gzz(i,j-1,k) ) * odysq180

      d2_gd(1,1,3,3) = (  2*gxx(i,j,k+3) - 27*gxx(i,j,k+2) + 270*gxx(i,j,k+1) - 490*gxx(i,j,k)&
                        + 2*gxx(i,j,k-3) - 27*gxx(i,j,k-2) + 270*gxx(i,j,k-1) ) * odzsq180
      d2_gd(1,2,3,3) = (  2*gxy(i,j,k+3) - 27*gxy(i,j,k+2) + 270*gxy(i,j,k+1) - 490*gxy(i,j,k)&
                        + 2*gxy(i,j,k-3) - 27*gxy(i,j,k-2) + 270*gxy(i,j,k-1) ) * odzsq180
      d2_gd(1,3,3,3) = (  2*gxz(i,j,k+3) - 27*gxz(i,j,k+2) + 270*gxz(i,j,k+1) - 490*gxz(i,j,k)&
                        + 2*gxz(i,j,k-3) - 27*gxz(i,j,k-2) + 270*gxz(i,j,k-1) ) * odzsq180
      d2_gd(2,2,3,3) = (  2*gyy(i,j,k+3) - 27*gyy(i,j,k+2) + 270*gyy(i,j,k+1) - 490*gyy(i,j,k)&
                        + 2*gyy(i,j,k-3) - 27*gyy(i,j,k-2) + 270*gyy(i,j,k-1) ) * odzsq180
      d2_gd(2,3,3,3) = (  2*gyz(i,j,k+3) - 27*gyz(i,j,k+2) + 270*gyz(i,j,k+1) - 490*gyz(i,j,k)&
                        + 2*gyz(i,j,k-3) - 27*gyz(i,j,k-2) + 270*gyz(i,j,k-1) ) * odzsq180
      d2_gd(3,3,3,3) = (  2*gzz(i,j,k+3) - 27*gzz(i,j,k+2) + 270*gzz(i,j,k+1) - 490*gzz(i,j,k)&
                        + 2*gzz(i,j,k-3) - 27*gzz(i,j,k-2) + 270*gzz(i,j,k-1) ) * odzsq180

      d2_gd(1,1,1,2) = (    -gxx(i-3,j+3,k) +   9*gxx(i-2,j+3,k) -   45*gxx(i-1,j+3,k) +   45*gxx(i+1,j+3,k) -   9*gxx(i+2,j+3,k) +&
  &    gxx(i+3,j+3,k) &
                        +  9*gxx(i-3,j+2,k) -  81*gxx(i-2,j+2,k) +  405*gxx(i-1,j+2,k) -  405*gxx(i+1,j+2,k) +  81*gxx(i+2,j+2,k) -&
  &  9*gxx(i+3,j+2,k) &
                        - 45*gxx(i-3,j+1,k) + 405*gxx(i-2,j+1,k) - 2025*gxx(i-1,j+1,k) + 2025*gxx(i+1,j+1,k) - 405*gxx(i+2,j+1,k) +&
  & 45*gxx(i+3,j+1,k) &
                        + 45*gxx(i-3,j-1,k) - 405*gxx(i-2,j-1,k) + 2025*gxx(i-1,j-1,k) - 2025*gxx(i+1,j-1,k) + 405*gxx(i+2,j-1,k) -&
  & 45*gxx(i+3,j-1,k) &
                        -  9*gxx(i-3,j-2,k) +  81*gxx(i-2,j-2,k) -  405*gxx(i-1,j-2,k) +  405*gxx(i+1,j-2,k) -  81*gxx(i+2,j-2,k) +&
  &  9*gxx(i+3,j-2,k) &
                        +    gxx(i-3,j-3,k) -   9*gxx(i-2,j-3,k) +   45*gxx(i-1,j-3,k) -   45*gxx(i+1,j-3,k) +   9*gxx(i+2,j-3,k) -&
  &    gxx(i+3,j-3,k) ) * odxdy3600
      d2_gd(1,2,1,2) = (    -gxy(i-3,j+3,k) +   9*gxy(i-2,j+3,k) -   45*gxy(i-1,j+3,k) +   45*gxy(i+1,j+3,k) -   9*gxy(i+2,j+3,k) +&
  &    gxy(i+3,j+3,k) &
                        +  9*gxy(i-3,j+2,k) -  81*gxy(i-2,j+2,k) +  405*gxy(i-1,j+2,k) -  405*gxy(i+1,j+2,k) +  81*gxy(i+2,j+2,k) -&
  &  9*gxy(i+3,j+2,k) &
                        - 45*gxy(i-3,j+1,k) + 405*gxy(i-2,j+1,k) - 2025*gxy(i-1,j+1,k) + 2025*gxy(i+1,j+1,k) - 405*gxy(i+2,j+1,k) +&
  & 45*gxy(i+3,j+1,k) &
                        + 45*gxy(i-3,j-1,k) - 405*gxy(i-2,j-1,k) + 2025*gxy(i-1,j-1,k) - 2025*gxy(i+1,j-1,k) + 405*gxy(i+2,j-1,k) -&
  & 45*gxy(i+3,j-1,k) &
                        -  9*gxy(i-3,j-2,k) +  81*gxy(i-2,j-2,k) -  405*gxy(i-1,j-2,k) +  405*gxy(i+1,j-2,k) -  81*gxy(i+2,j-2,k) +&
  &  9*gxy(i+3,j-2,k) &
                        +    gxy(i-3,j-3,k) -   9*gxy(i-2,j-3,k) +   45*gxy(i-1,j-3,k) -   45*gxy(i+1,j-3,k) +   9*gxy(i+2,j-3,k) -&
  &    gxy(i+3,j-3,k) ) * odxdy3600
      d2_gd(1,3,1,2) = (    -gxz(i-3,j+3,k) +   9*gxz(i-2,j+3,k) -   45*gxz(i-1,j+3,k) +   45*gxz(i+1,j+3,k) -   9*gxz(i+2,j+3,k) +&
  &    gxz(i+3,j+3,k) &
                        +  9*gxz(i-3,j+2,k) -  81*gxz(i-2,j+2,k) +  405*gxz(i-1,j+2,k) -  405*gxz(i+1,j+2,k) +  81*gxz(i+2,j+2,k) -&
  &  9*gxz(i+3,j+2,k) &
                        - 45*gxz(i-3,j+1,k) + 405*gxz(i-2,j+1,k) - 2025*gxz(i-1,j+1,k) + 2025*gxz(i+1,j+1,k) - 405*gxz(i+2,j+1,k) +&
  & 45*gxz(i+3,j+1,k) &
                        + 45*gxz(i-3,j-1,k) - 405*gxz(i-2,j-1,k) + 2025*gxz(i-1,j-1,k) - 2025*gxz(i+1,j-1,k) + 405*gxz(i+2,j-1,k) -&
  & 45*gxz(i+3,j-1,k) &
                        -  9*gxz(i-3,j-2,k) +  81*gxz(i-2,j-2,k) -  405*gxz(i-1,j-2,k) +  405*gxz(i+1,j-2,k) -  81*gxz(i+2,j-2,k) +&
  &  9*gxz(i+3,j-2,k) &
                        +    gxz(i-3,j-3,k) -   9*gxz(i-2,j-3,k) +   45*gxz(i-1,j-3,k) -   45*gxz(i+1,j-3,k) +   9*gxz(i+2,j-3,k) -&
  &    gxz(i+3,j-3,k) ) * odxdy3600
      d2_gd(2,2,1,2) = (    -gyy(i-3,j+3,k) +   9*gyy(i-2,j+3,k) -   45*gyy(i-1,j+3,k) +   45*gyy(i+1,j+3,k) -   9*gyy(i+2,j+3,k) +&
  &    gyy(i+3,j+3,k) &
                        +  9*gyy(i-3,j+2,k) -  81*gyy(i-2,j+2,k) +  405*gyy(i-1,j+2,k) -  405*gyy(i+1,j+2,k) +  81*gyy(i+2,j+2,k) -&
  &  9*gyy(i+3,j+2,k) &
                        - 45*gyy(i-3,j+1,k) + 405*gyy(i-2,j+1,k) - 2025*gyy(i-1,j+1,k) + 2025*gyy(i+1,j+1,k) - 405*gyy(i+2,j+1,k) +&
  & 45*gyy(i+3,j+1,k) &
                        + 45*gyy(i-3,j-1,k) - 405*gyy(i-2,j-1,k) + 2025*gyy(i-1,j-1,k) - 2025*gyy(i+1,j-1,k) + 405*gyy(i+2,j-1,k) -&
  & 45*gyy(i+3,j-1,k) &
                        -  9*gyy(i-3,j-2,k) +  81*gyy(i-2,j-2,k) -  405*gyy(i-1,j-2,k) +  405*gyy(i+1,j-2,k) -  81*gyy(i+2,j-2,k) +&
  &  9*gyy(i+3,j-2,k) &
                        +    gyy(i-3,j-3,k) -   9*gyy(i-2,j-3,k) +   45*gyy(i-1,j-3,k) -   45*gyy(i+1,j-3,k) +   9*gyy(i+2,j-3,k) -&
  &    gyy(i+3,j-3,k) ) * odxdy3600
      d2_gd(2,3,1,2) = (    -gyz(i-3,j+3,k) +   9*gyz(i-2,j+3,k) -   45*gyz(i-1,j+3,k) +   45*gyz(i+1,j+3,k) -   9*gyz(i+2,j+3,k) +&
  &    gyz(i+3,j+3,k) &
                        +  9*gyz(i-3,j+2,k) -  81*gyz(i-2,j+2,k) +  405*gyz(i-1,j+2,k) -  405*gyz(i+1,j+2,k) +  81*gyz(i+2,j+2,k) -&
  &  9*gyz(i+3,j+2,k) &
                        - 45*gyz(i-3,j+1,k) + 405*gyz(i-2,j+1,k) - 2025*gyz(i-1,j+1,k) + 2025*gyz(i+1,j+1,k) - 405*gyz(i+2,j+1,k) +&
  & 45*gyz(i+3,j+1,k) &
                        + 45*gyz(i-3,j-1,k) - 405*gyz(i-2,j-1,k) + 2025*gyz(i-1,j-1,k) - 2025*gyz(i+1,j-1,k) + 405*gyz(i+2,j-1,k) -&
  & 45*gyz(i+3,j-1,k) &
                        -  9*gyz(i-3,j-2,k) +  81*gyz(i-2,j-2,k) -  405*gyz(i-1,j-2,k) +  405*gyz(i+1,j-2,k) -  81*gyz(i+2,j-2,k) +&
  &  9*gyz(i+3,j-2,k) &
                        +    gyz(i-3,j-3,k) -   9*gyz(i-2,j-3,k) +   45*gyz(i-1,j-3,k) -   45*gyz(i+1,j-3,k) +   9*gyz(i+2,j-3,k) -&
  &    gyz(i+3,j-3,k) ) * odxdy3600
      d2_gd(3,3,1,2) = (    -gzz(i-3,j+3,k) +   9*gzz(i-2,j+3,k) -   45*gzz(i-1,j+3,k) +   45*gzz(i+1,j+3,k) -   9*gzz(i+2,j+3,k) +&
  &    gzz(i+3,j+3,k) &
                        +  9*gzz(i-3,j+2,k) -  81*gzz(i-2,j+2,k) +  405*gzz(i-1,j+2,k) -  405*gzz(i+1,j+2,k) +  81*gzz(i+2,j+2,k) -&
  &  9*gzz(i+3,j+2,k) &
                        - 45*gzz(i-3,j+1,k) + 405*gzz(i-2,j+1,k) - 2025*gzz(i-1,j+1,k) + 2025*gzz(i+1,j+1,k) - 405*gzz(i+2,j+1,k) +&
  & 45*gzz(i+3,j+1,k) &
                        + 45*gzz(i-3,j-1,k) - 405*gzz(i-2,j-1,k) + 2025*gzz(i-1,j-1,k) - 2025*gzz(i+1,j-1,k) + 405*gzz(i+2,j-1,k) -&
  & 45*gzz(i+3,j-1,k) &
                        -  9*gzz(i-3,j-2,k) +  81*gzz(i-2,j-2,k) -  405*gzz(i-1,j-2,k) +  405*gzz(i+1,j-2,k) -  81*gzz(i+2,j-2,k) +&
  &  9*gzz(i+3,j-2,k) &
                        +    gzz(i-3,j-3,k) -   9*gzz(i-2,j-3,k) +   45*gzz(i-1,j-3,k) -   45*gzz(i+1,j-3,k) +   9*gzz(i+2,j-3,k) -&
  &    gzz(i+3,j-3,k) ) * odxdy3600

      d2_gd(1,1,1,3) = (    -gxx(i-3,j,k+3) +   9*gxx(i-2,j,k+3) -   45*gxx(i-1,j,k+3) +   45*gxx(i+1,j,k+3) -   9*gxx(i+2,j,k+3) +&
  &    gxx(i+3,j,k+3) &
                        +  9*gxx(i-3,j,k+2) -  81*gxx(i-2,j,k+2) +  405*gxx(i-1,j,k+2) -  405*gxx(i+1,j,k+2) +  81*gxx(i+2,j,k+2) -&
  &  9*gxx(i+3,j,k+2) &
                        - 45*gxx(i-3,j,k+1) + 405*gxx(i-2,j,k+1) - 2025*gxx(i-1,j,k+1) + 2025*gxx(i+1,j,k+1) - 405*gxx(i+2,j,k+1) +&
  & 45*gxx(i+3,j,k+1) &
                        + 45*gxx(i-3,j,k-1) - 405*gxx(i-2,j,k-1) + 2025*gxx(i-1,j,k-1) - 2025*gxx(i+1,j,k-1) + 405*gxx(i+2,j,k-1) -&
  & 45*gxx(i+3,j,k-1) &
                        -  9*gxx(i-3,j,k-2) +  81*gxx(i-2,j,k-2) -  405*gxx(i-1,j,k-2) +  405*gxx(i+1,j,k-2) -  81*gxx(i+2,j,k-2) +&
  &  9*gxx(i+3,j,k-2) &
                        +    gxx(i-3,j,k-3) -   9*gxx(i-2,j,k-3) +   45*gxx(i-1,j,k-3) -   45*gxx(i+1,j,k-3) +   9*gxx(i+2,j,k-3) -&
  &    gxx(i+3,j,k-3) ) * odxdz3600
      d2_gd(1,2,1,3) = (    -gxy(i-3,j,k+3) +   9*gxy(i-2,j,k+3) -   45*gxy(i-1,j,k+3) +   45*gxy(i+1,j,k+3) -   9*gxy(i+2,j,k+3) +&
  &    gxy(i+3,j,k+3) &
                        +  9*gxy(i-3,j,k+2) -  81*gxy(i-2,j,k+2) +  405*gxy(i-1,j,k+2) -  405*gxy(i+1,j,k+2) +  81*gxy(i+2,j,k+2) -&
  &  9*gxy(i+3,j,k+2) &
                          - 45*gxy(i-3,j,k+1) + 405*gxy(i-2,j,k+1) - 2025*gxy(i-1,j,k+1) + 2025*gxy(i+1,j,k+1) - 405*gxy(i+2,j,k+1)&
  & + 45*gxy(i+3,j,k+1) &
                        + 45*gxy(i-3,j,k-1) - 405*gxy(i-2,j,k-1) + 2025*gxy(i-1,j,k-1) - 2025*gxy(i+1,j,k-1) + 405*gxy(i+2,j,k-1) -&
  & 45*gxy(i+3,j,k-1) &
                        -  9*gxy(i-3,j,k-2) +  81*gxy(i-2,j,k-2) -  405*gxy(i-1,j,k-2) +  405*gxy(i+1,j,k-2) -  81*gxy(i+2,j,k-2) +&
  &  9*gxy(i+3,j,k-2) &
                        +    gxy(i-3,j,k-3) -   9*gxy(i-2,j,k-3) +   45*gxy(i-1,j,k-3) -   45*gxy(i+1,j,k-3) +   9*gxy(i+2,j,k-3) -&
  &    gxy(i+3,j,k-3) ) * odxdz3600
      d2_gd(1,3,1,3) = (    -gxz(i-3,j,k+3) +   9*gxz(i-2,j,k+3) -   45*gxz(i-1,j,k+3) +   45*gxz(i+1,j,k+3) -   9*gxz(i+2,j,k+3) +&
  &    gxz(i+3,j,k+3) &
                        +  9*gxz(i-3,j,k+2) -  81*gxz(i-2,j,k+2) +  405*gxz(i-1,j,k+2) -  405*gxz(i+1,j,k+2) +  81*gxz(i+2,j,k+2) -&
  &  9*gxz(i+3,j,k+2) &
                        - 45*gxz(i-3,j,k+1) + 405*gxz(i-2,j,k+1) - 2025*gxz(i-1,j,k+1) + 2025*gxz(i+1,j,k+1) - 405*gxz(i+2,j,k+1) +&
  & 45*gxz(i+3,j,k+1) &
                        + 45*gxz(i-3,j,k-1) - 405*gxz(i-2,j,k-1) + 2025*gxz(i-1,j,k-1) - 2025*gxz(i+1,j,k-1) + 405*gxz(i+2,j,k-1) -&
  & 45*gxz(i+3,j,k-1) &
                        -  9*gxz(i-3,j,k-2) +  81*gxz(i-2,j,k-2) -  405*gxz(i-1,j,k-2) +  405*gxz(i+1,j,k-2) -  81*gxz(i+2,j,k-2) +&
  &  9*gxz(i+3,j,k-2) &
                        +    gxz(i-3,j,k-3) -   9*gxz(i-2,j,k-3) +   45*gxz(i-1,j,k-3) -   45*gxz(i+1,j,k-3) +   9*gxz(i+2,j,k-3) -&
  &    gxz(i+3,j,k-3) ) * odxdz3600
      d2_gd(2,2,1,3) = (    -gyy(i-3,j,k+3) +   9*gyy(i-2,j,k+3) -   45*gyy(i-1,j,k+3) +   45*gyy(i+1,j,k+3) -   9*gyy(i+2,j,k+3) +&
  &    gyy(i+3,j,k+3) &
                        +  9*gyy(i-3,j,k+2) -  81*gyy(i-2,j,k+2) +  405*gyy(i-1,j,k+2) -  405*gyy(i+1,j,k+2) +  81*gyy(i+2,j,k+2) -&
  &  9*gyy(i+3,j,k+2) &
                        - 45*gyy(i-3,j,k+1) + 405*gyy(i-2,j,k+1) - 2025*gyy(i-1,j,k+1) + 2025*gyy(i+1,j,k+1) - 405*gyy(i+2,j,k+1) +&
  & 45*gyy(i+3,j,k+1) &
                        + 45*gyy(i-3,j,k-1) - 405*gyy(i-2,j,k-1) + 2025*gyy(i-1,j,k-1) - 2025*gyy(i+1,j,k-1) + 405*gyy(i+2,j,k-1) -&
  & 45*gyy(i+3,j,k-1) &
                        -  9*gyy(i-3,j,k-2) +  81*gyy(i-2,j,k-2) -  405*gyy(i-1,j,k-2) +  405*gyy(i+1,j,k-2) -  81*gyy(i+2,j,k-2) +&
  &  9*gyy(i+3,j,k-2) &
                        +    gyy(i-3,j,k-3) -   9*gyy(i-2,j,k-3) +   45*gyy(i-1,j,k-3) -   45*gyy(i+1,j,k-3) +   9*gyy(i+2,j,k-3) -&
  &    gyy(i+3,j,k-3) ) * odxdz3600
      d2_gd(2,3,1,3) = (    -gyz(i-3,j,k+3) +   9*gyz(i-2,j,k+3) -   45*gyz(i-1,j,k+3) +   45*gyz(i+1,j,k+3) -   9*gyz(i+2,j,k+3) +&
  &    gyz(i+3,j,k+3) &
                        +  9*gyz(i-3,j,k+2) -  81*gyz(i-2,j,k+2) +  405*gyz(i-1,j,k+2) -  405*gyz(i+1,j,k+2) +  81*gyz(i+2,j,k+2) -&
  &  9*gyz(i+3,j,k+2) &
                        - 45*gyz(i-3,j,k+1) + 405*gyz(i-2,j,k+1) - 2025*gyz(i-1,j,k+1) + 2025*gyz(i+1,j,k+1) - 405*gyz(i+2,j,k+1) +&
  & 45*gyz(i+3,j,k+1) &
                        + 45*gyz(i-3,j,k-1) - 405*gyz(i-2,j,k-1) + 2025*gyz(i-1,j,k-1) - 2025*gyz(i+1,j,k-1) + 405*gyz(i+2,j,k-1) -&
  & 45*gyz(i+3,j,k-1) &
                        -  9*gyz(i-3,j,k-2) +  81*gyz(i-2,j,k-2) -  405*gyz(i-1,j,k-2) +  405*gyz(i+1,j,k-2) -  81*gyz(i+2,j,k-2) +&
  &  9*gyz(i+3,j,k-2) &
                        +    gyz(i-3,j,k-3) -   9*gyz(i-2,j,k-3) +   45*gyz(i-1,j,k-3) -   45*gyz(i+1,j,k-3) +   9*gyz(i+2,j,k-3) -&
  &    gyz(i+3,j,k-3) ) * odxdz3600
      d2_gd(3,3,1,3) = (    -gzz(i-3,j,k+3) +   9*gzz(i-2,j,k+3) -   45*gzz(i-1,j,k+3) +   45*gzz(i+1,j,k+3) -   9*gzz(i+2,j,k+3) +&
  &    gzz(i+3,j,k+3) &
                        +  9*gzz(i-3,j,k+2) -  81*gzz(i-2,j,k+2) +  405*gzz(i-1,j,k+2) -  405*gzz(i+1,j,k+2) +  81*gzz(i+2,j,k+2) -&
  &  9*gzz(i+3,j,k+2) &
                        - 45*gzz(i-3,j,k+1) + 405*gzz(i-2,j,k+1) - 2025*gzz(i-1,j,k+1) + 2025*gzz(i+1,j,k+1) - 405*gzz(i+2,j,k+1) +&
  & 45*gzz(i+3,j,k+1) &
                        + 45*gzz(i-3,j,k-1) - 405*gzz(i-2,j,k-1) + 2025*gzz(i-1,j,k-1) - 2025*gzz(i+1,j,k-1) + 405*gzz(i+2,j,k-1) -&
  & 45*gzz(i+3,j,k-1) &
                        -  9*gzz(i-3,j,k-2) +  81*gzz(i-2,j,k-2) -  405*gzz(i-1,j,k-2) +  405*gzz(i+1,j,k-2) -  81*gzz(i+2,j,k-2) +&
  &  9*gzz(i+3,j,k-2) &
                        +    gzz(i-3,j,k-3) -   9*gzz(i-2,j,k-3) +   45*gzz(i-1,j,k-3) -   45*gzz(i+1,j,k-3) +   9*gzz(i+2,j,k-3) -&
  &    gzz(i+3,j,k-3) ) * odxdz3600

      d2_gd(1,1,2,3) = (    -gxx(i,j-3,k+3) +   9*gxx(i,j-2,k+3) -   45*gxx(i,j-1,k+3) +   45*gxx(i,j+1,k+3) -   9*gxx(i,j+2,k+3) +&
  &    gxx(i,j+3,k+3) &
                        +  9*gxx(i,j-3,k+2) -  81*gxx(i,j-2,k+2) +  405*gxx(i,j-1,k+2) -  405*gxx(i,j+1,k+2) +  81*gxx(i,j+2,k+2) -&
  &  9*gxx(i,j+3,k+2) &
                        - 45*gxx(i,j-3,k+1) + 405*gxx(i,j-2,k+1) - 2025*gxx(i,j-1,k+1) + 2025*gxx(i,j+1,k+1) - 405*gxx(i,j+2,k+1) +&
  & 45*gxx(i,j+3,k+1) &
                        + 45*gxx(i,j-3,k-1) - 405*gxx(i,j-2,k-1) + 2025*gxx(i,j-1,k-1) - 2025*gxx(i,j+1,k-1) + 405*gxx(i,j+2,k-1) -&
  & 45*gxx(i,j+3,k-1) &
                        -  9*gxx(i,j-3,k-2) +  81*gxx(i,j-2,k-2) -  405*gxx(i,j-1,k-2) +  405*gxx(i,j+1,k-2) -  81*gxx(i,j+2,k-2) +&
  &  9*gxx(i,j+3,k-2) &
                        +    gxx(i,j-3,k-3) -   9*gxx(i,j-2,k-3) +   45*gxx(i,j-1,k-3) -   45*gxx(i,j+1,k-3) +   9*gxx(i,j+2,k-3) -&
  &    gxx(i,j+3,k-3) ) * odydz3600
      d2_gd(1,2,2,3) = (    -gxy(i,j-3,k+3) +   9*gxy(i,j-2,k+3) -   45*gxy(i,j-1,k+3) +   45*gxy(i,j+1,k+3) -   9*gxy(i,j+2,k+3) +&
  &    gxy(i,j+3,k+3) &
                        +  9*gxy(i,j-3,k+2) -  81*gxy(i,j-2,k+2) +  405*gxy(i,j-1,k+2) -  405*gxy(i,j+1,k+2) +  81*gxy(i,j+2,k+2) -&
  &  9*gxy(i,j+3,k+2) &
                        - 45*gxy(i,j-3,k+1) + 405*gxy(i,j-2,k+1) - 2025*gxy(i,j-1,k+1) + 2025*gxy(i,j+1,k+1) - 405*gxy(i,j+2,k+1) +&
  & 45*gxy(i,j+3,k+1) &
                        + 45*gxy(i,j-3,k-1) - 405*gxy(i,j-2,k-1) + 2025*gxy(i,j-1,k-1) - 2025*gxy(i,j+1,k-1) + 405*gxy(i,j+2,k-1) -&
  & 45*gxy(i,j+3,k-1) &
                        -  9*gxy(i,j-3,k-2) +  81*gxy(i,j-2,k-2) -  405*gxy(i,j-1,k-2) +  405*gxy(i,j+1,k-2) -  81*gxy(i,j+2,k-2) +&
  &  9*gxy(i,j+3,k-2) &
                        +    gxy(i,j-3,k-3) -   9*gxy(i,j-2,k-3) +   45*gxy(i,j-1,k-3) -   45*gxy(i,j+1,k-3) +   9*gxy(i,j+2,k-3) -&
  &    gxy(i,j+3,k-3) ) * odydz3600
      d2_gd(1,3,2,3) = (    -gxz(i,j-3,k+3) +   9*gxz(i,j-2,k+3) -   45*gxz(i,j-1,k+3) +   45*gxz(i,j+1,k+3) -   9*gxz(i,j+2,k+3) +&
  &    gxz(i,j+3,k+3) &
                        +  9*gxz(i,j-3,k+2) -  81*gxz(i,j-2,k+2) +  405*gxz(i,j-1,k+2) -  405*gxz(i,j+1,k+2) +  81*gxz(i,j+2,k+2) -&
  &  9*gxz(i,j+3,k+2) &
                        - 45*gxz(i,j-3,k+1) + 405*gxz(i,j-2,k+1) - 2025*gxz(i,j-1,k+1) + 2025*gxz(i,j+1,k+1) - 405*gxz(i,j+2,k+1) +&
  & 45*gxz(i,j+3,k+1) &
                        + 45*gxz(i,j-3,k-1) - 405*gxz(i,j-2,k-1) + 2025*gxz(i,j-1,k-1) - 2025*gxz(i,j+1,k-1) + 405*gxz(i,j+2,k-1) -&
  & 45*gxz(i,j+3,k-1) &
                        -  9*gxz(i,j-3,k-2) +  81*gxz(i,j-2,k-2) -  405*gxz(i,j-1,k-2) +  405*gxz(i,j+1,k-2) -  81*gxz(i,j+2,k-2) +&
  &  9*gxz(i,j+3,k-2) &
                        +    gxz(i,j-3,k-3) -   9*gxz(i,j-2,k-3) +   45*gxz(i,j-1,k-3) -   45*gxz(i,j+1,k-3) +   9*gxz(i,j+2,k-3) -&
  &    gxz(i,j+3,k-3) ) * odydz3600
      d2_gd(2,2,2,3) = (    -gyy(i,j-3,k+3) +   9*gyy(i,j-2,k+3) -   45*gyy(i,j-1,k+3) +   45*gyy(i,j+1,k+3) -   9*gyy(i,j+2,k+3) +&
  &    gyy(i,j+3,k+3) &
                        +  9*gyy(i,j-3,k+2) -  81*gyy(i,j-2,k+2) +  405*gyy(i,j-1,k+2) -  405*gyy(i,j+1,k+2) +  81*gyy(i,j+2,k+2) -&
  &  9*gyy(i,j+3,k+2) &
                        - 45*gyy(i,j-3,k+1) + 405*gyy(i,j-2,k+1) - 2025*gyy(i,j-1,k+1) + 2025*gyy(i,j+1,k+1) - 405*gyy(i,j+2,k+1) +&
  & 45*gyy(i,j+3,k+1) &
                        + 45*gyy(i,j-3,k-1) - 405*gyy(i,j-2,k-1) + 2025*gyy(i,j-1,k-1) - 2025*gyy(i,j+1,k-1) + 405*gyy(i,j+2,k-1) -&
  & 45*gyy(i,j+3,k-1) &
                        -  9*gyy(i,j-3,k-2) +  81*gyy(i,j-2,k-2) -  405*gyy(i,j-1,k-2) +  405*gyy(i,j+1,k-2) -  81*gyy(i,j+2,k-2) +&
  &  9*gyy(i,j+3,k-2) &
                        +    gyy(i,j-3,k-3) -   9*gyy(i,j-2,k-3) +   45*gyy(i,j-1,k-3) -   45*gyy(i,j+1,k-3) +   9*gyy(i,j+2,k-3) -&
  &    gyy(i,j+3,k-3) ) * odydz3600
      d2_gd(2,3,2,3) = (    -gyz(i,j-3,k+3) +   9*gyz(i,j-2,k+3) -   45*gyz(i,j-1,k+3) +   45*gyz(i,j+1,k+3) -   9*gyz(i,j+2,k+3) +&
  &    gyz(i,j+3,k+3) &
                        +  9*gyz(i,j-3,k+2) -  81*gyz(i,j-2,k+2) +  405*gyz(i,j-1,k+2) -  405*gyz(i,j+1,k+2) +  81*gyz(i,j+2,k+2) -&
  &  9*gyz(i,j+3,k+2) &
                        - 45*gyz(i,j-3,k+1) + 405*gyz(i,j-2,k+1) - 2025*gyz(i,j-1,k+1) + 2025*gyz(i,j+1,k+1) - 405*gyz(i,j+2,k+1) +&
  & 45*gyz(i,j+3,k+1) &
                        + 45*gyz(i,j-3,k-1) - 405*gyz(i,j-2,k-1) + 2025*gyz(i,j-1,k-1) - 2025*gyz(i,j+1,k-1) + 405*gyz(i,j+2,k-1) -&
  & 45*gyz(i,j+3,k-1) &
                        -  9*gyz(i,j-3,k-2) +  81*gyz(i,j-2,k-2) -  405*gyz(i,j-1,k-2) +  405*gyz(i,j+1,k-2) -  81*gyz(i,j+2,k-2) +&
  &  9*gyz(i,j+3,k-2) &
                        +    gyz(i,j-3,k-3) -   9*gyz(i,j-2,k-3) +   45*gyz(i,j-1,k-3) -   45*gyz(i,j+1,k-3) +   9*gyz(i,j+2,k-3) -&
  &    gyz(i,j+3,k-3) ) * odydz3600
      d2_gd(3,3,2,3) = (    -gzz(i,j-3,k+3) +   9*gzz(i,j-2,k+3) -   45*gzz(i,j-1,k+3) +   45*gzz(i,j+1,k+3) -   9*gzz(i,j+2,k+3) +&
  &    gzz(i,j+3,k+3) &
                        +  9*gzz(i,j-3,k+2) -  81*gzz(i,j-2,k+2) +  405*gzz(i,j-1,k+2) -  405*gzz(i,j+1,k+2) +  81*gzz(i,j+2,k+2) -&
  &  9*gzz(i,j+3,k+2) &
                        - 45*gzz(i,j-3,k+1) + 405*gzz(i,j-2,k+1) - 2025*gzz(i,j-1,k+1) + 2025*gzz(i,j+1,k+1) - 405*gzz(i,j+2,k+1) +&
  & 45*gzz(i,j+3,k+1) &
                        + 45*gzz(i,j-3,k-1) - 405*gzz(i,j-2,k-1) + 2025*gzz(i,j-1,k-1) - 2025*gzz(i,j+1,k-1) + 405*gzz(i,j+2,k-1) -&
  & 45*gzz(i,j+3,k-1) &
                        -  9*gzz(i,j-3,k-2) +  81*gzz(i,j-2,k-2) -  405*gzz(i,j-1,k-2) +  405*gzz(i,j+1,k-2) -  81*gzz(i,j+2,k-2) +&
  &  9*gzz(i,j+3,k-2) &
                        +    gzz(i,j-3,k-3) -   9*gzz(i,j-2,k-3) +   45*gzz(i,j-1,k-3) -   45*gzz(i,j+1,k-3) +   9*gzz(i,j+2,k-3) -&
  &    gzz(i,j+3,k-3) ) * odydz3600

      d2_gd(1,1,2,1) = d2_gd(1,1,1,2)
      d2_gd(1,2,2,1) = d2_gd(1,2,1,2)
      d2_gd(1,3,2,1) = d2_gd(1,3,1,2)
      d2_gd(2,2,2,1) = d2_gd(2,2,1,2)
      d2_gd(2,3,2,1) = d2_gd(2,3,1,2)
      d2_gd(3,3,2,1) = d2_gd(3,3,1,2)

      d2_gd(1,1,3,1) = d2_gd(1,1,1,3)
      d2_gd(1,2,3,1) = d2_gd(1,2,1,3)
      d2_gd(1,3,3,1) = d2_gd(1,3,1,3)
      d2_gd(2,2,3,1) = d2_gd(2,2,1,3)
      d2_gd(2,3,3,1) = d2_gd(2,3,1,3)
      d2_gd(3,3,3,1) = d2_gd(3,3,1,3)

      d2_gd(1,1,3,2) = d2_gd(1,1,2,3)
      d2_gd(1,2,3,2) = d2_gd(1,2,2,3)
      d2_gd(1,3,3,2) = d2_gd(1,3,2,3)
      d2_gd(2,2,3,2) = d2_gd(2,2,2,3)
      d2_gd(2,3,3,2) = d2_gd(2,3,2,3)
      d2_gd(3,3,3,2) = d2_gd(3,3,2,3)

      d2_gd(2,1,:,:) = d2_gd(1,2,:,:)
      d2_gd(3,1,:,:) = d2_gd(1,3,:,:)
      d2_gd(3,2,:,:) = d2_gd(2,3,:,:)


    
    else if ( NP_order == 4 ) then
    
    
    
    d1_gd(1,1,1) = (   -gxx(i+2,j,k) + 8*gxx(i+1,j,k)                      &
                    - 8*gxx(i-1,j,k) +   gxx(i-2,j,k) ) / dx12
    d1_gd(1,2,1) = (   -gxy(i+2,j,k) + 8*gxy(i+1,j,k)                      &
                    - 8*gxy(i-1,j,k) +   gxy(i-2,j,k) ) / dx12
    d1_gd(1,3,1) = (   -gxz(i+2,j,k) + 8*gxz(i+1,j,k)                      &
                    - 8*gxz(i-1,j,k) +   gxz(i-2,j,k) ) / dx12
    d1_gd(2,2,1) = (   -gyy(i+2,j,k) + 8*gyy(i+1,j,k)                      &
                    - 8*gyy(i-1,j,k) +   gyy(i-2,j,k) ) / dx12
    d1_gd(2,3,1) = (   -gyz(i+2,j,k) + 8*gyz(i+1,j,k)                      &
                    - 8*gyz(i-1,j,k) +   gyz(i-2,j,k) ) / dx12
    d1_gd(3,3,1) = (   -gzz(i+2,j,k) + 8*gzz(i+1,j,k)                      &
                    - 8*gzz(i-1,j,k) +   gzz(i-2,j,k) ) / dx12

    d1_gd(1,1,2) = (   -gxx(i,j+2,k) + 8*gxx(i,j+1,k)                      &
                    - 8*gxx(i,j-1,k) +   gxx(i,j-2,k) ) / dy12
    d1_gd(1,2,2) = (   -gxy(i,j+2,k) + 8*gxy(i,j+1,k)                      &
                    - 8*gxy(i,j-1,k) +   gxy(i,j-2,k) ) / dy12
    d1_gd(1,3,2) = (   -gxz(i,j+2,k) + 8*gxz(i,j+1,k)                      &
                    - 8*gxz(i,j-1,k) +   gxz(i,j-2,k) ) / dy12
    d1_gd(2,2,2) = (   -gyy(i,j+2,k) + 8*gyy(i,j+1,k)                      &
                    - 8*gyy(i,j-1,k) +   gyy(i,j-2,k) ) / dy12
    d1_gd(2,3,2) = (   -gyz(i,j+2,k) + 8*gyz(i,j+1,k)                      &
                    - 8*gyz(i,j-1,k) +   gyz(i,j-2,k) ) / dy12
    d1_gd(3,3,2) = (   -gzz(i,j+2,k) + 8*gzz(i,j+1,k)                      &
                    - 8*gzz(i,j-1,k) +   gzz(i,j-2,k) ) / dy12

    d1_gd(1,1,3) = (   -gxx(i,j,k+2) + 8*gxx(i,j,k+1)                      &
                    - 8*gxx(i,j,k-1) +   gxx(i,j,k-2) ) / dz12
    d1_gd(1,2,3) = (   -gxy(i,j,k+2) + 8*gxy(i,j,k+1)                      &
                    - 8*gxy(i,j,k-1) +   gxy(i,j,k-2) ) / dz12
    d1_gd(1,3,3) = (   -gxz(i,j,k+2) + 8*gxz(i,j,k+1)                      &
                    - 8*gxz(i,j,k-1) +   gxz(i,j,k-2) ) / dz12
    d1_gd(2,2,3) = (   -gyy(i,j,k+2) + 8*gyy(i,j,k+1)                      &
                    - 8*gyy(i,j,k-1) +   gyy(i,j,k-2) ) / dz12
    d1_gd(2,3,3) = (   -gyz(i,j,k+2) + 8*gyz(i,j,k+1)                      &
                    - 8*gyz(i,j,k-1) +   gyz(i,j,k-2) ) / dz12
    d1_gd(3,3,3) = (   -gzz(i,j,k+2) + 8*gzz(i,j,k+1)                      &
                    - 8*gzz(i,j,k-1) +   gzz(i,j,k-2) ) / dz12

    d1_gd(2,1,:) = d1_gd(1,2,:)
    d1_gd(3,1,:) = d1_gd(1,3,:)
    d1_gd(3,2,:) = d1_gd(2,3,:)


    
    d1_kd(1,1,1) = (   -kxx(i+2,j,k) + 8*kxx(i+1,j,k)                      &
                    - 8*kxx(i-1,j,k) +   kxx(i-2,j,k) ) / dx12
    d1_kd(1,2,1) = (   -kxy(i+2,j,k) + 8*kxy(i+1,j,k)                      &
                    - 8*kxy(i-1,j,k) +   kxy(i-2,j,k) ) / dx12
    d1_kd(1,3,1) = (   -kxz(i+2,j,k) + 8*kxz(i+1,j,k)                      &
                    - 8*kxz(i-1,j,k) +   kxz(i-2,j,k) ) / dx12
    d1_kd(2,2,1) = (   -kyy(i+2,j,k) + 8*kyy(i+1,j,k)                      &
                    - 8*kyy(i-1,j,k) +   kyy(i-2,j,k) ) / dx12
    d1_kd(2,3,1) = (   -kyz(i+2,j,k) + 8*kyz(i+1,j,k)                      &
                    - 8*kyz(i-1,j,k) +   kyz(i-2,j,k) ) / dx12
    d1_kd(3,3,1) = (   -kzz(i+2,j,k) + 8*kzz(i+1,j,k)                      &
                    - 8*kzz(i-1,j,k) +   kzz(i-2,j,k) ) / dx12

    d1_kd(1,1,2) = (   -kxx(i,j+2,k) + 8*kxx(i,j+1,k)                      &
                    - 8*kxx(i,j-1,k) +   kxx(i,j-2,k) ) / dy12
    d1_kd(1,2,2) = (   -kxy(i,j+2,k) + 8*kxy(i,j+1,k)                      &
                    - 8*kxy(i,j-1,k) +   kxy(i,j-2,k) ) / dy12
    d1_kd(1,3,2) = (   -kxz(i,j+2,k) + 8*kxz(i,j+1,k)                      &
                    - 8*kxz(i,j-1,k) +   kxz(i,j-2,k) ) / dy12
    d1_kd(2,2,2) = (   -kyy(i,j+2,k) + 8*kyy(i,j+1,k)                      &
                    - 8*kyy(i,j-1,k) +   kyy(i,j-2,k) ) / dy12
    d1_kd(2,3,2) = (   -kyz(i,j+2,k) + 8*kyz(i,j+1,k)                      &
                    - 8*kyz(i,j-1,k) +   kyz(i,j-2,k) ) / dy12
    d1_kd(3,3,2) = (   -kzz(i,j+2,k) + 8*kzz(i,j+1,k)                      &
                    - 8*kzz(i,j-1,k) +   kzz(i,j-2,k) ) / dy12

    d1_kd(1,1,3) = (   -kxx(i,j,k+2) + 8*kxx(i,j,k+1)                      &
                    - 8*kxx(i,j,k-1) +   kxx(i,j,k-2) ) / dz12
    d1_kd(1,2,3) = (   -kxy(i,j,k+2) + 8*kxy(i,j,k+1)                      &
                    - 8*kxy(i,j,k-1) +   kxy(i,j,k-2) ) / dz12
    d1_kd(1,3,3) = (   -kxz(i,j,k+2) + 8*kxz(i,j,k+1)                      &
                    - 8*kxz(i,j,k-1) +   kxz(i,j,k-2) ) / dz12
    d1_kd(2,2,3) = (   -kyy(i,j,k+2) + 8*kyy(i,j,k+1)                      &
                    - 8*kyy(i,j,k-1) +   kyy(i,j,k-2) ) / dz12
    d1_kd(2,3,3) = (   -kyz(i,j,k+2) + 8*kyz(i,j,k+1)                      &
                    - 8*kyz(i,j,k-1) +   kyz(i,j,k-2) ) / dz12
    d1_kd(3,3,3) = (   -kzz(i,j,k+2) + 8*kzz(i,j,k+1)                      &
                    - 8*kzz(i,j,k-1) +   kzz(i,j,k-2) ) / dz12

    d1_kd(2,1,:) = d1_kd(1,2,:)
    d1_kd(3,1,:) = d1_kd(1,3,:)
    d1_kd(3,2,:) = d1_kd(2,3,:)
    

    
    
    d2_gd(1,1,1,1) = (   -gxx(i+2,j,k) + 16*gxx(i+1,j,k) - 30*gxx(i,j,k)   &
                     + 16*gxx(i-1,j,k) -    gxx(i-2,j,k) ) / dxsq12
    d2_gd(1,2,1,1) = (   -gxy(i+2,j,k) + 16*gxy(i+1,j,k) - 30*gxy(i,j,k)   &
                     + 16*gxy(i-1,j,k) -    gxy(i-2,j,k) ) / dxsq12
    d2_gd(1,3,1,1) = (   -gxz(i+2,j,k) + 16*gxz(i+1,j,k) - 30*gxz(i,j,k)   &
                     + 16*gxz(i-1,j,k) -    gxz(i-2,j,k) ) / dxsq12
    d2_gd(2,2,1,1) = (   -gyy(i+2,j,k) + 16*gyy(i+1,j,k) - 30*gyy(i,j,k)   &
                     + 16*gyy(i-1,j,k) -    gyy(i-2,j,k) ) / dxsq12
    d2_gd(2,3,1,1) = (   -gyz(i+2,j,k) + 16*gyz(i+1,j,k) - 30*gyz(i,j,k)   &
                     + 16*gyz(i-1,j,k) -    gyz(i-2,j,k) ) / dxsq12
    d2_gd(3,3,1,1) = (   -gzz(i+2,j,k) + 16*gzz(i+1,j,k) - 30*gzz(i,j,k)   &
                     + 16*gzz(i-1,j,k) -    gzz(i-2,j,k) ) / dxsq12

    d2_gd(1,1,2,2) = (   -gxx(i,j+2,k) + 16*gxx(i,j+1,k) - 30*gxx(i,j,k)   &
                     + 16*gxx(i,j-1,k) -    gxx(i,j-2,k) ) / dysq12
    d2_gd(1,2,2,2) = (   -gxy(i,j+2,k) + 16*gxy(i,j+1,k) - 30*gxy(i,j,k)   &
                     + 16*gxy(i,j-1,k) -    gxy(i,j-2,k) ) / dysq12
    d2_gd(1,3,2,2) = (   -gxz(i,j+2,k) + 16*gxz(i,j+1,k) - 30*gxz(i,j,k)   &
                     + 16*gxz(i,j-1,k) -    gxz(i,j-2,k) ) / dysq12
    d2_gd(2,2,2,2) = (   -gyy(i,j+2,k) + 16*gyy(i,j+1,k) - 30*gyy(i,j,k)   &
                     + 16*gyy(i,j-1,k) -    gyy(i,j-2,k) ) / dysq12
    d2_gd(2,3,2,2) = (   -gyz(i,j+2,k) + 16*gyz(i,j+1,k) - 30*gyz(i,j,k)   &
                     + 16*gyz(i,j-1,k) -    gyz(i,j-2,k) ) / dysq12
    d2_gd(3,3,2,2) = (   -gzz(i,j+2,k) + 16*gzz(i,j+1,k) - 30*gzz(i,j,k)   &
                     + 16*gzz(i,j-1,k) -    gzz(i,j-2,k) ) / dysq12

    d2_gd(1,1,3,3) = (   -gxx(i,j,k+2) + 16*gxx(i,j,k+1) - 30*gxx(i,j,k)   &
                     + 16*gxx(i,j,k-1) -    gxx(i,j,k-2) ) / dzsq12
    d2_gd(1,2,3,3) = (   -gxy(i,j,k+2) + 16*gxy(i,j,k+1) - 30*gxy(i,j,k)   &
                     + 16*gxy(i,j,k-1) -    gxy(i,j,k-2) ) / dzsq12
    d2_gd(1,3,3,3) = (   -gxz(i,j,k+2) + 16*gxz(i,j,k+1) - 30*gxz(i,j,k)   &
                     + 16*gxz(i,j,k-1) -    gxz(i,j,k-2) ) / dzsq12
    d2_gd(2,2,3,3) = (   -gyy(i,j,k+2) + 16*gyy(i,j,k+1) - 30*gyy(i,j,k)   &
                     + 16*gyy(i,j,k-1) -    gyy(i,j,k-2) ) / dzsq12
    d2_gd(2,3,3,3) = (   -gyz(i,j,k+2) + 16*gyz(i,j,k+1) - 30*gyz(i,j,k)   &
                     + 16*gyz(i,j,k-1) -    gyz(i,j,k-2) ) / dzsq12
    d2_gd(3,3,3,3) = (   -gzz(i,j,k+2) + 16*gzz(i,j,k+1) - 30*gzz(i,j,k)   &
                     + 16*gzz(i,j,k-1) -    gzz(i,j,k-2) ) / dzsq12

    d2_gd(1,1,1,2) = (   -gxx(i-2,j+2,k) +  8*gxx(i-1,j+2,k) -  8*gxx(i+1,j+2,k) +   gxx(i+2,j+2,k)   &
                      + 8*gxx(i-2,j+1,k) - 64*gxx(i-1,j+1,k) + 64*gxx(i+1,j+1,k) - 8*gxx(i+2,j+1,k)   &
                      - 8*gxx(i-2,j-1,k) + 64*gxx(i-1,j-1,k) - 64*gxx(i+1,j-1,k) + 8*gxx(i+2,j-1,k)   &
                      +   gxx(i-2,j-2,k) -  8*gxx(i-1,j-2,k) +  8*gxx(i+1,j-2,k) -   gxx(i+2,j-2,k) ) / dxdy144
    d2_gd(1,2,1,2) = (   -gxy(i-2,j+2,k) +  8*gxy(i-1,j+2,k) -  8*gxy(i+1,j+2,k) +   gxy(i+2,j+2,k)   &
                      + 8*gxy(i-2,j+1,k) - 64*gxy(i-1,j+1,k) + 64*gxy(i+1,j+1,k) - 8*gxy(i+2,j+1,k)   &
                      - 8*gxy(i-2,j-1,k) + 64*gxy(i-1,j-1,k) - 64*gxy(i+1,j-1,k) + 8*gxy(i+2,j-1,k)   &
                      +   gxy(i-2,j-2,k) -  8*gxy(i-1,j-2,k) +  8*gxy(i+1,j-2,k) -   gxy(i+2,j-2,k) ) / dxdy144
    d2_gd(1,3,1,2) = (   -gxz(i-2,j+2,k) +  8*gxz(i-1,j+2,k) -  8*gxz(i+1,j+2,k) +   gxz(i+2,j+2,k)   &
                      + 8*gxz(i-2,j+1,k) - 64*gxz(i-1,j+1,k) + 64*gxz(i+1,j+1,k) - 8*gxz(i+2,j+1,k)   &
                      - 8*gxz(i-2,j-1,k) + 64*gxz(i-1,j-1,k) - 64*gxz(i+1,j-1,k) + 8*gxz(i+2,j-1,k)   &
                      +   gxz(i-2,j-2,k) -  8*gxz(i-1,j-2,k) +  8*gxz(i+1,j-2,k) -   gxz(i+2,j-2,k) ) / dxdy144
    d2_gd(2,2,1,2) = (   -gyy(i-2,j+2,k) +  8*gyy(i-1,j+2,k) -  8*gyy(i+1,j+2,k) +   gyy(i+2,j+2,k)   &
                      + 8*gyy(i-2,j+1,k) - 64*gyy(i-1,j+1,k) + 64*gyy(i+1,j+1,k) - 8*gyy(i+2,j+1,k)   &
                      - 8*gyy(i-2,j-1,k) + 64*gyy(i-1,j-1,k) - 64*gyy(i+1,j-1,k) + 8*gyy(i+2,j-1,k)   &
                      +   gyy(i-2,j-2,k) -  8*gyy(i-1,j-2,k) +  8*gyy(i+1,j-2,k) -   gyy(i+2,j-2,k) ) / dxdy144
    d2_gd(2,3,1,2) = (   -gyz(i-2,j+2,k) +  8*gyz(i-1,j+2,k) -  8*gyz(i+1,j+2,k) +   gyz(i+2,j+2,k)   &
                      + 8*gyz(i-2,j+1,k) - 64*gyz(i-1,j+1,k) + 64*gyz(i+1,j+1,k) - 8*gyz(i+2,j+1,k)   &
                      - 8*gyz(i-2,j-1,k) + 64*gyz(i-1,j-1,k) - 64*gyz(i+1,j-1,k) + 8*gyz(i+2,j-1,k)   &
                      +   gyz(i-2,j-2,k) -  8*gyz(i-1,j-2,k) +  8*gyz(i+1,j-2,k) -   gyz(i+2,j-2,k) ) / dxdy144
    d2_gd(3,3,1,2) = (   -gzz(i-2,j+2,k) +  8*gzz(i-1,j+2,k) -  8*gzz(i+1,j+2,k) +   gzz(i+2,j+2,k)   &
                      + 8*gzz(i-2,j+1,k) - 64*gzz(i-1,j+1,k) + 64*gzz(i+1,j+1,k) - 8*gzz(i+2,j+1,k)   &
                      - 8*gzz(i-2,j-1,k) + 64*gzz(i-1,j-1,k) - 64*gzz(i+1,j-1,k) + 8*gzz(i+2,j-1,k)   &
                      +   gzz(i-2,j-2,k) -  8*gzz(i-1,j-2,k) +  8*gzz(i+1,j-2,k) -   gzz(i+2,j-2,k) ) / dxdy144

    d2_gd(1,1,1,3) = (   -gxx(i-2,j,k+2) +  8*gxx(i-1,j,k+2) -  8*gxx(i+1,j,k+2) +   gxx(i+2,j,k+2)   &
                      + 8*gxx(i-2,j,k+1) - 64*gxx(i-1,j,k+1) + 64*gxx(i+1,j,k+1) - 8*gxx(i+2,j,k+1)   &
                      - 8*gxx(i-2,j,k-1) + 64*gxx(i-1,j,k-1) - 64*gxx(i+1,j,k-1) + 8*gxx(i+2,j,k-1)   &
                      +   gxx(i-2,j,k-2) -  8*gxx(i-1,j,k-2) +  8*gxx(i+1,j,k-2) -   gxx(i+2,j,k-2) ) / dxdz144
    d2_gd(1,2,1,3) = (   -gxy(i-2,j,k+2) +  8*gxy(i-1,j,k+2) -  8*gxy(i+1,j,k+2) +   gxy(i+2,j,k+2)   &
                      + 8*gxy(i-2,j,k+1) - 64*gxy(i-1,j,k+1) + 64*gxy(i+1,j,k+1) - 8*gxy(i+2,j,k+1)   &
                      - 8*gxy(i-2,j,k-1) + 64*gxy(i-1,j,k-1) - 64*gxy(i+1,j,k-1) + 8*gxy(i+2,j,k-1)   &
                      +   gxy(i-2,j,k-2) -  8*gxy(i-1,j,k-2) +  8*gxy(i+1,j,k-2) -   gxy(i+2,j,k-2) ) / dxdz144
    d2_gd(1,3,1,3) = (   -gxz(i-2,j,k+2) +  8*gxz(i-1,j,k+2) -  8*gxz(i+1,j,k+2) +   gxz(i+2,j,k+2)   &
                      + 8*gxz(i-2,j,k+1) - 64*gxz(i-1,j,k+1) + 64*gxz(i+1,j,k+1) - 8*gxz(i+2,j,k+1)   &
                      - 8*gxz(i-2,j,k-1) + 64*gxz(i-1,j,k-1) - 64*gxz(i+1,j,k-1) + 8*gxz(i+2,j,k-1)   &
                      +   gxz(i-2,j,k-2) -  8*gxz(i-1,j,k-2) +  8*gxz(i+1,j,k-2) -   gxz(i+2,j,k-2) ) / dxdz144
    d2_gd(2,2,1,3) = (   -gyy(i-2,j,k+2) +  8*gyy(i-1,j,k+2) -  8*gyy(i+1,j,k+2) +   gyy(i+2,j,k+2)   &
                      + 8*gyy(i-2,j,k+1) - 64*gyy(i-1,j,k+1) + 64*gyy(i+1,j,k+1) - 8*gyy(i+2,j,k+1)   &
                      - 8*gyy(i-2,j,k-1) + 64*gyy(i-1,j,k-1) - 64*gyy(i+1,j,k-1) + 8*gyy(i+2,j,k-1)   &
                      +   gyy(i-2,j,k-2) -  8*gyy(i-1,j,k-2) +  8*gyy(i+1,j,k-2) -   gyy(i+2,j,k-2) ) / dxdz144
    d2_gd(2,3,1,3) = (   -gyz(i-2,j,k+2) +  8*gyz(i-1,j,k+2) -  8*gyz(i+1,j,k+2) +   gyz(i+2,j,k+2)   &
                      + 8*gyz(i-2,j,k+1) - 64*gyz(i-1,j,k+1) + 64*gyz(i+1,j,k+1) - 8*gyz(i+2,j,k+1)   &
                      - 8*gyz(i-2,j,k-1) + 64*gyz(i-1,j,k-1) - 64*gyz(i+1,j,k-1) + 8*gyz(i+2,j,k-1)   &
                      +   gyz(i-2,j,k-2) -  8*gyz(i-1,j,k-2) +  8*gyz(i+1,j,k-2) -   gyz(i+2,j,k-2) ) / dxdz144
    d2_gd(3,3,1,3) = (   -gzz(i-2,j,k+2) +  8*gzz(i-1,j,k+2) -  8*gzz(i+1,j,k+2) +   gzz(i+2,j,k+2)   &
                      + 8*gzz(i-2,j,k+1) - 64*gzz(i-1,j,k+1) + 64*gzz(i+1,j,k+1) - 8*gzz(i+2,j,k+1)   &
                      - 8*gzz(i-2,j,k-1) + 64*gzz(i-1,j,k-1) - 64*gzz(i+1,j,k-1) + 8*gzz(i+2,j,k-1)   &
                      +   gzz(i-2,j,k-2) -  8*gzz(i-1,j,k-2) +  8*gzz(i+1,j,k-2) -   gzz(i+2,j,k-2) ) / dxdz144

    d2_gd(1,1,2,3) = (   -gxx(i,j-2,k+2) +  8*gxx(i,j-1,k+2) -  8*gxx(i,j+1,k+2) +   gxx(i,j+2,k+2)   &
                      + 8*gxx(i,j-2,k+1) - 64*gxx(i,j-1,k+1) + 64*gxx(i,j+1,k+1) - 8*gxx(i,j+2,k+1)   &
                      - 8*gxx(i,j-2,k-1) + 64*gxx(i,j-1,k-1) - 64*gxx(i,j+1,k-1) + 8*gxx(i,j+2,k-1)   &
                      +   gxx(i,j-2,k-2) -  8*gxx(i,j-1,k-2) +  8*gxx(i,j+1,k-2) -   gxx(i,j+2,k-2) ) / dydz144
    d2_gd(1,2,2,3) = (   -gxy(i,j-2,k+2) +  8*gxy(i,j-1,k+2) -  8*gxy(i,j+1,k+2) +   gxy(i,j+2,k+2)   &
                      + 8*gxy(i,j-2,k+1) - 64*gxy(i,j-1,k+1) + 64*gxy(i,j+1,k+1) - 8*gxy(i,j+2,k+1)   &
                      - 8*gxy(i,j-2,k-1) + 64*gxy(i,j-1,k-1) - 64*gxy(i,j+1,k-1) + 8*gxy(i,j+2,k-1)   &
                      +   gxy(i,j-2,k-2) -  8*gxy(i,j-1,k-2) +  8*gxy(i,j+1,k-2) -   gxy(i,j+2,k-2) ) / dydz144
    d2_gd(1,3,2,3) = (   -gxz(i,j-2,k+2) +  8*gxz(i,j-1,k+2) -  8*gxz(i,j+1,k+2) +   gxz(i,j+2,k+2)   &
                      + 8*gxz(i,j-2,k+1) - 64*gxz(i,j-1,k+1) + 64*gxz(i,j+1,k+1) - 8*gxz(i,j+2,k+1)   &
                      - 8*gxz(i,j-2,k-1) + 64*gxz(i,j-1,k-1) - 64*gxz(i,j+1,k-1) + 8*gxz(i,j+2,k-1)   &
                      +   gxz(i,j-2,k-2) -  8*gxz(i,j-1,k-2) +  8*gxz(i,j+1,k-2) -   gxz(i,j+2,k-2) ) / dydz144
    d2_gd(2,2,2,3) = (   -gyy(i,j-2,k+2) +  8*gyy(i,j-1,k+2) -  8*gyy(i,j+1,k+2) +   gyy(i,j+2,k+2)   &
                      + 8*gyy(i,j-2,k+1) - 64*gyy(i,j-1,k+1) + 64*gyy(i,j+1,k+1) - 8*gyy(i,j+2,k+1)   &
                      - 8*gyy(i,j-2,k-1) + 64*gyy(i,j-1,k-1) - 64*gyy(i,j+1,k-1) + 8*gyy(i,j+2,k-1)   &
                      +   gyy(i,j-2,k-2) -  8*gyy(i,j-1,k-2) +  8*gyy(i,j+1,k-2) -   gyy(i,j+2,k-2) ) / dydz144
    d2_gd(2,3,2,3) = (   -gyz(i,j-2,k+2) +  8*gyz(i,j-1,k+2) -  8*gyz(i,j+1,k+2) +   gyz(i,j+2,k+2)   &
                      + 8*gyz(i,j-2,k+1) - 64*gyz(i,j-1,k+1) + 64*gyz(i,j+1,k+1) - 8*gyz(i,j+2,k+1)   &
                      - 8*gyz(i,j-2,k-1) + 64*gyz(i,j-1,k-1) - 64*gyz(i,j+1,k-1) + 8*gyz(i,j+2,k-1)   &
                      +   gyz(i,j-2,k-2) -  8*gyz(i,j-1,k-2) +  8*gyz(i,j+1,k-2) -   gyz(i,j+2,k-2) ) / dydz144
    d2_gd(3,3,2,3) = (   -gzz(i,j-2,k+2) +  8*gzz(i,j-1,k+2) -  8*gzz(i,j+1,k+2) +   gzz(i,j+2,k+2)   &
                      + 8*gzz(i,j-2,k+1) - 64*gzz(i,j-1,k+1) + 64*gzz(i,j+1,k+1) - 8*gzz(i,j+2,k+1)   &
                      - 8*gzz(i,j-2,k-1) + 64*gzz(i,j-1,k-1) - 64*gzz(i,j+1,k-1) + 8*gzz(i,j+2,k-1)   &
                      +   gzz(i,j-2,k-2) -  8*gzz(i,j-1,k-2) +  8*gzz(i,j+1,k-2) -   gzz(i,j+2,k-2) ) / dydz144

    d2_gd(1,1,2,1) = d2_gd(1,1,1,2)
    d2_gd(1,2,2,1) = d2_gd(1,2,1,2)
    d2_gd(1,3,2,1) = d2_gd(1,3,1,2)
    d2_gd(2,2,2,1) = d2_gd(2,2,1,2)
    d2_gd(2,3,2,1) = d2_gd(2,3,1,2)
    d2_gd(3,3,2,1) = d2_gd(3,3,1,2)

    d2_gd(1,1,3,1) = d2_gd(1,1,1,3)
    d2_gd(1,2,3,1) = d2_gd(1,2,1,3)
    d2_gd(1,3,3,1) = d2_gd(1,3,1,3)
    d2_gd(2,2,3,1) = d2_gd(2,2,1,3)
    d2_gd(2,3,3,1) = d2_gd(2,3,1,3)
    d2_gd(3,3,3,1) = d2_gd(3,3,1,3)

    d2_gd(1,1,3,2) = d2_gd(1,1,2,3)
    d2_gd(1,2,3,2) = d2_gd(1,2,2,3)
    d2_gd(1,3,3,2) = d2_gd(1,3,2,3)
    d2_gd(2,2,3,2) = d2_gd(2,2,2,3)
    d2_gd(2,3,3,2) = d2_gd(2,3,2,3)
    d2_gd(3,3,3,2) = d2_gd(3,3,2,3)

    d2_gd(2,1,:,:) = d2_gd(1,2,:,:)
    d2_gd(3,1,:,:) = d2_gd(1,3,:,:)
    d2_gd(3,2,:,:) = d2_gd(2,3,:,:)
    

    else
    
    
    
    d1_gd(1,1,1) = (gxx(i+1,j,k) - gxx(i-1,j,k)) / dx2
    d1_gd(1,2,1) = (gxy(i+1,j,k) - gxy(i-1,j,k)) / dx2
    d1_gd(1,3,1) = (gxz(i+1,j,k) - gxz(i-1,j,k)) / dx2
    d1_gd(2,2,1) = (gyy(i+1,j,k) - gyy(i-1,j,k)) / dx2
    d1_gd(2,3,1) = (gyz(i+1,j,k) - gyz(i-1,j,k)) / dx2
    d1_gd(3,3,1) = (gzz(i+1,j,k) - gzz(i-1,j,k)) / dx2

    d1_gd(1,1,2) = (gxx(i,j+1,k) - gxx(i,j-1,k)) / dy2
    d1_gd(1,2,2) = (gxy(i,j+1,k) - gxy(i,j-1,k)) / dy2
    d1_gd(1,3,2) = (gxz(i,j+1,k) - gxz(i,j-1,k)) / dy2
    d1_gd(2,2,2) = (gyy(i,j+1,k) - gyy(i,j-1,k)) / dy2
    d1_gd(2,3,2) = (gyz(i,j+1,k) - gyz(i,j-1,k)) / dy2
    d1_gd(3,3,2) = (gzz(i,j+1,k) - gzz(i,j-1,k)) / dy2

    d1_gd(1,1,3) = (gxx(i,j,k+1) - gxx(i,j,k-1)) / dz2
    d1_gd(1,2,3) = (gxy(i,j,k+1) - gxy(i,j,k-1)) / dz2
    d1_gd(1,3,3) = (gxz(i,j,k+1) - gxz(i,j,k-1)) / dz2
    d1_gd(2,2,3) = (gyy(i,j,k+1) - gyy(i,j,k-1)) / dz2
    d1_gd(2,3,3) = (gyz(i,j,k+1) - gyz(i,j,k-1)) / dz2
    d1_gd(3,3,3) = (gzz(i,j,k+1) - gzz(i,j,k-1)) / dz2

    d1_gd(2,1,:) = d1_gd(1,2,:)
    d1_gd(3,1,:) = d1_gd(1,3,:)
    d1_gd(3,2,:) = d1_gd(2,3,:)


    
    d1_kd(1,1,1) = (kxx(i+1,j,k) - kxx(i-1,j,k)) / dx2
    d1_kd(1,2,1) = (kxy(i+1,j,k) - kxy(i-1,j,k)) / dx2
    d1_kd(1,3,1) = (kxz(i+1,j,k) - kxz(i-1,j,k)) / dx2
    d1_kd(2,2,1) = (kyy(i+1,j,k) - kyy(i-1,j,k)) / dx2
    d1_kd(2,3,1) = (kyz(i+1,j,k) - kyz(i-1,j,k)) / dx2
    d1_kd(3,3,1) = (kzz(i+1,j,k) - kzz(i-1,j,k)) / dx2

    d1_kd(1,1,2) = (kxx(i,j+1,k) - kxx(i,j-1,k)) / dy2
    d1_kd(1,2,2) = (kxy(i,j+1,k) - kxy(i,j-1,k)) / dy2
    d1_kd(1,3,2) = (kxz(i,j+1,k) - kxz(i,j-1,k)) / dy2
    d1_kd(2,2,2) = (kyy(i,j+1,k) - kyy(i,j-1,k)) / dy2
    d1_kd(2,3,2) = (kyz(i,j+1,k) - kyz(i,j-1,k)) / dy2
    d1_kd(3,3,2) = (kzz(i,j+1,k) - kzz(i,j-1,k)) / dy2

    d1_kd(1,1,3) = (kxx(i,j,k+1) - kxx(i,j,k-1)) / dz2
    d1_kd(1,2,3) = (kxy(i,j,k+1) - kxy(i,j,k-1)) / dz2
    d1_kd(1,3,3) = (kxz(i,j,k+1) - kxz(i,j,k-1)) / dz2
    d1_kd(2,2,3) = (kyy(i,j,k+1) - kyy(i,j,k-1)) / dz2
    d1_kd(2,3,3) = (kyz(i,j,k+1) - kyz(i,j,k-1)) / dz2
    d1_kd(3,3,3) = (kzz(i,j,k+1) - kzz(i,j,k-1)) / dz2

    d1_kd(2,1,:) = d1_kd(1,2,:)
    d1_kd(3,1,:) = d1_kd(1,3,:)
    d1_kd(3,2,:) = d1_kd(2,3,:)

    
    
    d2_gd(1,1,1,1) = (gxx(i+1,j,k) - 2*gxx(i,j,k) + gxx(i-1,j,k)) / dxsq
    d2_gd(1,2,1,1) = (gxy(i+1,j,k) - 2*gxy(i,j,k) + gxy(i-1,j,k)) / dxsq
    d2_gd(1,3,1,1) = (gxz(i+1,j,k) - 2*gxz(i,j,k) + gxz(i-1,j,k)) / dxsq
    d2_gd(2,2,1,1) = (gyy(i+1,j,k) - 2*gyy(i,j,k) + gyy(i-1,j,k)) / dxsq
    d2_gd(2,3,1,1) = (gyz(i+1,j,k) - 2*gyz(i,j,k) + gyz(i-1,j,k)) / dxsq
    d2_gd(3,3,1,1) = (gzz(i+1,j,k) - 2*gzz(i,j,k) + gzz(i-1,j,k)) / dxsq

    d2_gd(1,1,2,2) = (gxx(i,j+1,k) - 2*gxx(i,j,k) + gxx(i,j-1,k)) / dysq
    d2_gd(1,2,2,2) = (gxy(i,j+1,k) - 2*gxy(i,j,k) + gxy(i,j-1,k)) / dysq
    d2_gd(1,3,2,2) = (gxz(i,j+1,k) - 2*gxz(i,j,k) + gxz(i,j-1,k)) / dysq
    d2_gd(2,2,2,2) = (gyy(i,j+1,k) - 2*gyy(i,j,k) + gyy(i,j-1,k)) / dysq
    d2_gd(2,3,2,2) = (gyz(i,j+1,k) - 2*gyz(i,j,k) + gyz(i,j-1,k)) / dysq
    d2_gd(3,3,2,2) = (gzz(i,j+1,k) - 2*gzz(i,j,k) + gzz(i,j-1,k)) / dysq

    d2_gd(1,1,3,3) = (gxx(i,j,k+1) - 2*gxx(i,j,k) + gxx(i,j,k-1)) / dzsq
    d2_gd(1,2,3,3) = (gxy(i,j,k+1) - 2*gxy(i,j,k) + gxy(i,j,k-1)) / dzsq
    d2_gd(1,3,3,3) = (gxz(i,j,k+1) - 2*gxz(i,j,k) + gxz(i,j,k-1)) / dzsq
    d2_gd(2,2,3,3) = (gyy(i,j,k+1) - 2*gyy(i,j,k) + gyy(i,j,k-1)) / dzsq
    d2_gd(2,3,3,3) = (gyz(i,j,k+1) - 2*gyz(i,j,k) + gyz(i,j,k-1)) / dzsq
    d2_gd(3,3,3,3) = (gzz(i,j,k+1) - 2*gzz(i,j,k) + gzz(i,j,k-1)) / dzsq

    d2_gd(1,1,1,2) = ( gxx(i+1,j+1,k) + gxx(i-1,j-1,k)                     &
                     - gxx(i+1,j-1,k) - gxx(i-1,j+1,k) ) / dxdy4
    d2_gd(1,2,1,2) = ( gxy(i+1,j+1,k) + gxy(i-1,j-1,k)                     &
                     - gxy(i+1,j-1,k) - gxy(i-1,j+1,k) ) / dxdy4
    d2_gd(1,3,1,2) = ( gxz(i+1,j+1,k) + gxz(i-1,j-1,k)                     &
                     - gxz(i+1,j-1,k) - gxz(i-1,j+1,k) ) / dxdy4
    d2_gd(2,2,1,2) = ( gyy(i+1,j+1,k) + gyy(i-1,j-1,k)                     &
                     - gyy(i+1,j-1,k) - gyy(i-1,j+1,k) ) / dxdy4
    d2_gd(2,3,1,2) = ( gyz(i+1,j+1,k) + gyz(i-1,j-1,k)                     &
                     - gyz(i+1,j-1,k) - gyz(i-1,j+1,k) ) / dxdy4
    d2_gd(3,3,1,2) = ( gzz(i+1,j+1,k) + gzz(i-1,j-1,k)                     &
                     - gzz(i+1,j-1,k) - gzz(i-1,j+1,k) ) / dxdy4

    d2_gd(1,1,1,3) = ( gxx(i+1,j,k+1) + gxx(i-1,j,k-1)                     &
                     - gxx(i+1,j,k-1) - gxx(i-1,j,k+1) ) / dxdz4
    d2_gd(1,2,1,3) = ( gxy(i+1,j,k+1) + gxy(i-1,j,k-1)                     &
                     - gxy(i+1,j,k-1) - gxy(i-1,j,k+1) ) / dxdz4
    d2_gd(1,3,1,3) = ( gxz(i+1,j,k+1) + gxz(i-1,j,k-1)                     &
                     - gxz(i+1,j,k-1) - gxz(i-1,j,k+1) ) / dxdz4
    d2_gd(2,2,1,3) = ( gyy(i+1,j,k+1) + gyy(i-1,j,k-1)                     &
                     - gyy(i+1,j,k-1) - gyy(i-1,j,k+1) ) / dxdz4
    d2_gd(2,3,1,3) = ( gyz(i+1,j,k+1) + gyz(i-1,j,k-1)                     &
                     - gyz(i+1,j,k-1) - gyz(i-1,j,k+1) ) / dxdz4
    d2_gd(3,3,1,3) = ( gzz(i+1,j,k+1) + gzz(i-1,j,k-1)                     &
                     - gzz(i+1,j,k-1) - gzz(i-1,j,k+1) ) / dxdz4

    d2_gd(1,1,2,3) = ( gxx(i,j+1,k+1) + gxx(i,j-1,k-1)                     &
                     - gxx(i,j+1,k-1) - gxx(i,j-1,k+1) ) / dydz4
    d2_gd(1,2,2,3) = ( gxy(i,j+1,k+1) + gxy(i,j-1,k-1)                     &
                     - gxy(i,j+1,k-1) - gxy(i,j-1,k+1) ) / dydz4
    d2_gd(1,3,2,3) = ( gxz(i,j+1,k+1) + gxz(i,j-1,k-1)                     &
                     - gxz(i,j+1,k-1) - gxz(i,j-1,k+1) ) / dydz4
    d2_gd(2,2,2,3) = ( gyy(i,j+1,k+1) + gyy(i,j-1,k-1)                     &
                     - gyy(i,j+1,k-1) - gyy(i,j-1,k+1) ) / dydz4
    d2_gd(2,3,2,3) = ( gyz(i,j+1,k+1) + gyz(i,j-1,k-1)                     &
                     - gyz(i,j+1,k-1) - gyz(i,j-1,k+1) ) / dydz4
    d2_gd(3,3,2,3) = ( gzz(i,j+1,k+1) + gzz(i,j-1,k-1)                     &
                     - gzz(i,j+1,k-1) - gzz(i,j-1,k+1) ) / dydz4

    d2_gd(1,1,2,1) = d2_gd(1,1,1,2)
    d2_gd(1,2,2,1) = d2_gd(1,2,1,2)
    d2_gd(1,3,2,1) = d2_gd(1,3,1,2)
    d2_gd(2,2,2,1) = d2_gd(2,2,1,2)
    d2_gd(2,3,2,1) = d2_gd(2,3,1,2)
    d2_gd(3,3,2,1) = d2_gd(3,3,1,2)

    d2_gd(1,1,3,1) = d2_gd(1,1,1,3)
    d2_gd(1,2,3,1) = d2_gd(1,2,1,3)
    d2_gd(1,3,3,1) = d2_gd(1,3,1,3)
    d2_gd(2,2,3,1) = d2_gd(2,2,1,3)
    d2_gd(2,3,3,1) = d2_gd(2,3,1,3)
    d2_gd(3,3,3,1) = d2_gd(3,3,1,3)

    d2_gd(1,1,3,2) = d2_gd(1,1,2,3)
    d2_gd(1,2,3,2) = d2_gd(1,2,2,3)
    d2_gd(1,3,3,2) = d2_gd(1,3,2,3)
    d2_gd(2,2,3,2) = d2_gd(2,2,2,3)
    d2_gd(2,3,3,2) = d2_gd(2,3,2,3)
    d2_gd(3,3,3,2) = d2_gd(3,3,2,3)

    d2_gd(2,1,:,:) = d2_gd(1,2,:,:)
    d2_gd(3,1,:,:) = d2_gd(1,3,:,:)
    d2_gd(3,2,:,:) = d2_gd(2,3,:,:)

    end if
    

    
    cf1 = 0
    do a = 1, 3
      do b = 1, 3
        do c = 1, 3
          cf1(a,b,c) = 0.5d0 * (d1_gd(a,b,c) + d1_gd(a,c,b) - d1_gd(b,c,a))
        end do
      end do
    end do

    cf2 = 0
    do a = 1, 3
      do b = 1, 3
        do c = 1, 3
          do m = 1, 3
            cf2(a,b,c) = cf2(a,b,c) + gu(a,m) * cf1(m,b,c)
          end do
        end do
      end do
    end do
    


    
    cd1_kd = d1_kd
    do a = 1, 3
      do b = 1, 3
        do c = 1, 3
          do m = 1, 3
            cd1_kd(a,b,c) = cd1_kd(a,b,c) - cf2(m,a,c) * kd(m,b)           &
                                          - cf2(m,b,c) * kd(a,m)
          end do
        end do
      end do
    end do
    


    
    d1_gu = 0
    do a = 1, 3
      do b = 1, 3
        do c = 1, 3
          do m = 1, 3
            do n = 1, 3
              d1_gu(a,b,c) = d1_gu(a,b,c) - gu(a,m) * gu(b,n) * d1_gd(m,n,c)
            end do
          end do
        end do
      end do
    end do

    d1_cf2 = 0
    do a = 1, 3
      do b = 1, 3
        do c = 1, 3
          do d = 1, 3
            do m = 1, 3
              d1_cf2(a,b,c,d) = d1_cf2(a,b,c,d) + gu(a,m) * (d2_gd(c,m,b,d)&
                                + d2_gd(m,b,c,d) - d2_gd(b,c,m,d)) / 2     &
                                + d1_gu(a,m,d) * cf1(m,b,c)
            end do
          end do
        end do
      end do
    end do
    


    
    ri = 0
    do a = 1, 3
      do b = 1, 3
        do m = 1, 3
          ri(a,b) = ri(a,b) + d1_cf2(m,b,a,m) - d1_cf2(m,m,a,b)
          do n = 1, 3
            ri(a,b) = ri(a,b) + cf2(m,m,n) * cf2(n,b,a)                    &
                              - cf2(m,b,n) * cf2(n,m,a)
          end do
        end do
      end do
    end do
    


    
    eps_lc_u        = 0
    eps_lc_u(1,2,3) = 1
    eps_lc_u(2,3,1) = 1
    eps_lc_u(3,1,2) = 1
    eps_lc_u(3,2,1) = -1
    eps_lc_u(2,1,3) = -1
    eps_lc_u(1,3,2) = -1
    eps_lc_u = eps_lc_u / sqrt(detgd)

    eps_lc_d = eps_lc_u * detgd
    


    
    
    xx(:) = (/ x(i,j,k), y(i,j,k), z(i,j,k) /)

    
    
    

    if( xx(1)**2 + xx(2)**2 < 1.0d-12 ) xx(1) = xx(1) + 1.0d-10
    u_vec(:) = xx(:)
    v_vec(:) = (/ xx(1)*xx(3), xx(2)*xx(3), -xx(1)**2 - xx(2)**2 /)
    w_vec(:) = (/ -xx(2), xx(1), 0.0d0 /)

    
    dotp1 =   gd(1,1) * u_vec(1) * u_vec(1) + gd(1,2) * u_vec(1) * u_vec(2)&
            + gd(1,3) * u_vec(1) * u_vec(3) + gd(2,1) * u_vec(2) * u_vec(1)&
            + gd(2,2) * u_vec(2) * u_vec(2) + gd(2,3) * u_vec(2) * u_vec(3)&
            + gd(3,1) * u_vec(3) * u_vec(1) + gd(3,2) * u_vec(3) * u_vec(2)&
            + gd(3,3) * u_vec(3) * u_vec(3)
    u_vec = u_vec / sqrt(dotp1)

    dotp1 =   gd(1,1) * u_vec(1) * v_vec(1) + gd(1,2) * u_vec(1) * v_vec(2)&
            + gd(1,3) * u_vec(1) * v_vec(3) + gd(2,1) * u_vec(2) * v_vec(1)&
            + gd(2,2) * u_vec(2) * v_vec(2) + gd(2,3) * u_vec(2) * v_vec(3)&
            + gd(3,1) * u_vec(3) * v_vec(1) + gd(3,2) * u_vec(3) * v_vec(2)&
            + gd(3,3) * u_vec(3) * v_vec(3)
    v_vec = v_vec - dotp1 * u_vec

    dotp1 =   gd(1,1) * v_vec(1) * v_vec(1) + gd(1,2) * v_vec(1) * v_vec(2)&
            + gd(1,3) * v_vec(1) * v_vec(3) + gd(2,1) * v_vec(2) * v_vec(1)&
            + gd(2,2) * v_vec(2) * v_vec(2) + gd(2,3) * v_vec(2) * v_vec(3)&
            + gd(3,1) * v_vec(3) * v_vec(1) + gd(3,2) * v_vec(3) * v_vec(2)&
            + gd(3,3) * v_vec(3) * v_vec(3)
    v_vec = v_vec / sqrt(dotp1)

    dotp1 =   gd(1,1) * u_vec(1) * w_vec(1) + gd(1,2) * u_vec(1) * w_vec(2)&
            + gd(1,3) * u_vec(1) * w_vec(3) + gd(2,1) * u_vec(2) * w_vec(1)&
            + gd(2,2) * u_vec(2) * w_vec(2) + gd(2,3) * u_vec(2) * w_vec(3)&
            + gd(3,1) * u_vec(3) * w_vec(1) + gd(3,2) * u_vec(3) * w_vec(2)&
            + gd(3,3) * u_vec(3) * w_vec(3)

    dotp2 =   gd(1,1) * v_vec(1) * w_vec(1) + gd(1,2) * v_vec(1) * w_vec(2)&
            + gd(1,3) * v_vec(1) * w_vec(3) + gd(2,1) * v_vec(2) * w_vec(1)&
            + gd(2,2) * v_vec(2) * w_vec(2) + gd(2,3) * v_vec(2) * w_vec(3)&
            + gd(3,1) * v_vec(3) * w_vec(1) + gd(3,2) * v_vec(3) * w_vec(2)&
            + gd(3,3) * v_vec(3) * w_vec(3)
    w_vec = w_vec - dotp1 * u_vec - dotp2 * v_vec

    dotp1 =   gd(1,1) * w_vec(1) * w_vec(1) + gd(1,2) * w_vec(1) * w_vec(2)&
            + gd(1,3) * w_vec(1) * w_vec(3) + gd(2,1) * w_vec(2) * w_vec(1)&
            + gd(2,2) * w_vec(2) * w_vec(2) + gd(2,3) * w_vec(2) * w_vec(3)&
            + gd(3,1) * w_vec(3) * w_vec(1) + gd(3,2) * w_vec(3) * w_vec(2)&
            + gd(3,3) * w_vec(3) * w_vec(3)
    w_vec = w_vec / sqrt(dotp1)

    ud_vec = matmul( gd, u_vec )
    


    
    elec = ri
    mag  = 0
    do a = 1, 3
      do b = 1, 3
        do m = 1, 3
          do n = 1, 3
            elec(a,b) = elec(a,b) + gu(m,n) * (kd(a,b) * kd(m,n)           &
                                               - kd(a,m) * kd(b,n))
            do p = 1, 3
              mag(a,b) = mag(a,b) + gd(b,p) * eps_lc_u(p,m,n) * cd1_kd(n,a,m)   &
                                  + gd(a,p) * eps_lc_u(p,m,n) * cd1_kd(n,b,m)
            end do
          end do
        end do
      end do
    end do
    mag = 0.5 * mag

    
    electr = 0.0
    magtr  = 0.0
    do m = 1, 3
      do n = 1, 3
        electr = electr + gu(m,n) * elec(m,n)
        magtr  = magtr  + gu(m,n) * mag(m,n)
      end do
    end do

    elec = elec - gd * electr / 3.0
    mag  = mag  - gd * magtr  / 3.0

    


    
    psi4re(i,j,k) = 0
    psi4im(i,j,k) = 0
    do m = 1, 3
      do n = 1, 3
        psi4re(i,j,k) = psi4re(i,j,k)                                      &
                        + 0.5 * (   elec(m,n) * (   v_vec(m) * v_vec(n)    &
                                                  - w_vec(m) * w_vec(n) )  &
                                  - mag(m,n)  * (   v_vec(m) * w_vec(n)    &
                                                  + w_vec(m) * v_vec(n) ) )
        psi4im(i,j,k) = psi4im(i,j,k)                                      &
                        + 0.5 * (  -elec(m,n) * (   v_vec(m) * w_vec(n)    &
                                                  + w_vec(m) * v_vec(n) )  &
                                  + mag(m,n)  * (   w_vec(m) * w_vec(n)    &
                                                  - v_vec(m) * v_vec(n) ) )
      end do
    end do
    


  end do
  end do
  end do

end subroutine NP_calcPsiGF
