
























































































































 


































































































































































 



























































































































































































































 



















































 










 








































































































































































































































 










                                
                                

                                
                                

                                
                                

                                



















































































subroutine LeanBSSN_bssn_constraints( cctk_dim,cctk_gsh,cctk_lsh,cctk_lbnd,cctk_ubnd,cctk_tile_min,cctk_tile_max,cctk_ash,cctk_alig&
  &nment,cctk_alignment_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_ash&
  &2,cctk_ash3, active_slicing_handle,alp,alp_p,alp_p_p,axx,axx_p,axx_p_p,axy,axy_p,axy_p_p,axz,axz_p,axz_p_p,ayy,ayy_p,ayy_p_p,ayz&
  &,ayz_p,ayz_p_p,azz,azz_p,azz_p_p,betax,betax_p,betax_p_p,betay,betay_p,betay_p_p,betaz,betaz_p,betaz_p_p,coarse_dx,coarse_dy,coa&
  &rse_dz,conf_fac,conf_fac_p,conf_fac_p_p,conformal_state,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,eTtt,eTtt_p,eTtt_p_p,eTtx,eTtx_p,eTtx_p_p,eTty,eTty_p,eTty&
  &_p_p,eTtz,eTtz_p,eTtz_p_p,eTxx,eTxx_p,eTxx_p_p,eTxy,eTxy_p,eTxy_p_p,eTxz,eTxz_p,eTxz_p_p,eTyy,eTyy_p,eTyy_p_p,eTyz,eTyz_p,eTyz_p&
  &_p,eTzz,eTzz_p,eTzz_p_p,gammatx,gammatx_p,gammatx_p_p,gammaty,gammaty_p,gammaty_p_p,gammatz,gammatz_p,gammatz_p_p,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,hc,hc_p,hc_p_p,hxx,hxx_p,hxx_p_p,h&
  &xy,hxy_p,hxy_p_p,hxz,hxz_p,hxz_p_p,hyy,hyy_p,hyy_p_p,hyz,hyz_p,hyz_p_p,hzz,hzz_p,hzz_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,mcx,mcx_p,mcx_p_p,mcy,mcy_p,mcy_p_p,mcz,mcz_p,mcz_p_p,psi,p&
  &six,psixx,psixy,psixz,psiy,psiyy,psiyz,psiz,psizz,r,rhs_alp,rhs_axx,rhs_axy,rhs_axz,rhs_ayy,rhs_ayz,rhs_azz,rhs_betax,rhs_betay,&
  &rhs_betaz,rhs_conf_fac,rhs_gammatx,rhs_gammatx_p,rhs_gammatx_p_p,rhs_gammaty,rhs_gammaty_p,rhs_gammaty_p_p,rhs_gammatz,rhs_gamma&
  &tz_p,rhs_gammatz_p_p,rhs_hxx,rhs_hxy,rhs_hxz,rhs_hyy,rhs_hyz,rhs_hzz,rhs_tracek,shift_state,stress_energy_state,tracek,tracek_p,&
  &tracek_p_p,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)
  INTEGER*4 active_slicing_handle
  integer, parameter :: cctki_use_active_slicing_handle = kind(active_slicing_handle)
  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 axx (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_axx = kind(axx)
  REAL*8 axx_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_axx_p = kind(axx_p)
  REAL*8 axx_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_axx_p_p = kind(axx_p_p)
  REAL*8 axy (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_axy = kind(axy)
  REAL*8 axy_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_axy_p = kind(axy_p)
  REAL*8 axy_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_axy_p_p = kind(axy_p_p)
  REAL*8 axz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_axz = kind(axz)
  REAL*8 axz_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_axz_p = kind(axz_p)
  REAL*8 axz_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_axz_p_p = kind(axz_p_p)
  REAL*8 ayy (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_ayy = kind(ayy)
  REAL*8 ayy_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_ayy_p = kind(ayy_p)
  REAL*8 ayy_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_ayy_p_p = kind(ayy_p_p)
  REAL*8 ayz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_ayz = kind(ayz)
  REAL*8 ayz_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_ayz_p = kind(ayz_p)
  REAL*8 ayz_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_ayz_p_p = kind(ayz_p_p)
  REAL*8 azz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_azz = kind(azz)
  REAL*8 azz_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_azz_p = kind(azz_p)
  REAL*8 azz_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_azz_p_p = kind(azz_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 conf_fac (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_conf_fac = kind(conf_fac)
  REAL*8 conf_fac_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_conf_fac_p = kind(conf_fac_p)
  REAL*8 conf_fac_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_conf_fac_p_p = kind(conf_fac_p_p)
  INTEGER*4 conformal_state
  integer, parameter :: cctki_use_conformal_state = kind(conformal_state)
  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 eTtt (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_eTtt = kind(eTtt)
  REAL*8 eTtt_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_eTtt_p = kind(eTtt_p)
  REAL*8 eTtt_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_eTtt_p_p = kind(eTtt_p_p)
  REAL*8 eTtx (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_eTtx = kind(eTtx)
  REAL*8 eTtx_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_eTtx_p = kind(eTtx_p)
  REAL*8 eTtx_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_eTtx_p_p = kind(eTtx_p_p)
  REAL*8 eTty (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_eTty = kind(eTty)
  REAL*8 eTty_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_eTty_p = kind(eTty_p)
  REAL*8 eTty_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_eTty_p_p = kind(eTty_p_p)
  REAL*8 eTtz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_eTtz = kind(eTtz)
  REAL*8 eTtz_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_eTtz_p = kind(eTtz_p)
  REAL*8 eTtz_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_eTtz_p_p = kind(eTtz_p_p)
  REAL*8 eTxx (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_eTxx = kind(eTxx)
  REAL*8 eTxx_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_eTxx_p = kind(eTxx_p)
  REAL*8 eTxx_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_eTxx_p_p = kind(eTxx_p_p)
  REAL*8 eTxy (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_eTxy = kind(eTxy)
  REAL*8 eTxy_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_eTxy_p = kind(eTxy_p)
  REAL*8 eTxy_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_eTxy_p_p = kind(eTxy_p_p)
  REAL*8 eTxz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_eTxz = kind(eTxz)
  REAL*8 eTxz_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_eTxz_p = kind(eTxz_p)
  REAL*8 eTxz_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_eTxz_p_p = kind(eTxz_p_p)
  REAL*8 eTyy (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_eTyy = kind(eTyy)
  REAL*8 eTyy_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_eTyy_p = kind(eTyy_p)
  REAL*8 eTyy_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_eTyy_p_p = kind(eTyy_p_p)
  REAL*8 eTyz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_eTyz = kind(eTyz)
  REAL*8 eTyz_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_eTyz_p = kind(eTyz_p)
  REAL*8 eTyz_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_eTyz_p_p = kind(eTyz_p_p)
  REAL*8 eTzz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_eTzz = kind(eTzz)
  REAL*8 eTzz_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_eTzz_p = kind(eTzz_p)
  REAL*8 eTzz_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_eTzz_p_p = kind(eTzz_p_p)
  REAL*8 gammatx (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_gammatx = kind(gammatx)
  REAL*8 gammatx_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_gammatx_p = kind(gammatx_p)
  REAL*8 gammatx_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_gammatx_p_p = kind(gammatx_p_p)
  REAL*8 gammaty (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_gammaty = kind(gammaty)
  REAL*8 gammaty_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_gammaty_p = kind(gammaty_p)
  REAL*8 gammaty_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_gammaty_p_p = kind(gammaty_p_p)
  REAL*8 gammatz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_gammatz = kind(gammatz)
  REAL*8 gammatz_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_gammatz_p = kind(gammatz_p)
  REAL*8 gammatz_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_gammatz_p_p = kind(gammatz_p_p)
  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 hc (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_hc = kind(hc)
  REAL*8 hc_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_hc_p = kind(hc_p)
  REAL*8 hc_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_hc_p_p = kind(hc_p_p)
  REAL*8 hxx (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_hxx = kind(hxx)
  REAL*8 hxx_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_hxx_p = kind(hxx_p)
  REAL*8 hxx_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_hxx_p_p = kind(hxx_p_p)
  REAL*8 hxy (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_hxy = kind(hxy)
  REAL*8 hxy_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_hxy_p = kind(hxy_p)
  REAL*8 hxy_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_hxy_p_p = kind(hxy_p_p)
  REAL*8 hxz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_hxz = kind(hxz)
  REAL*8 hxz_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_hxz_p = kind(hxz_p)
  REAL*8 hxz_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_hxz_p_p = kind(hxz_p_p)
  REAL*8 hyy (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_hyy = kind(hyy)
  REAL*8 hyy_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_hyy_p = kind(hyy_p)
  REAL*8 hyy_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_hyy_p_p = kind(hyy_p_p)
  REAL*8 hyz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_hyz = kind(hyz)
  REAL*8 hyz_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_hyz_p = kind(hyz_p)
  REAL*8 hyz_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_hyz_p_p = kind(hyz_p_p)
  REAL*8 hzz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_hzz = kind(hzz)
  REAL*8 hzz_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_hzz_p = kind(hzz_p)
  REAL*8 hzz_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_hzz_p_p = kind(hzz_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 mcx (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_mcx = kind(mcx)
  REAL*8 mcx_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_mcx_p = kind(mcx_p)
  REAL*8 mcx_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_mcx_p_p = kind(mcx_p_p)
  REAL*8 mcy (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_mcy = kind(mcy)
  REAL*8 mcy_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_mcy_p = kind(mcy_p)
  REAL*8 mcy_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_mcy_p_p = kind(mcy_p_p)
  REAL*8 mcz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_mcz = kind(mcz)
  REAL*8 mcz_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_mcz_p = kind(mcz_p)
  REAL*8 mcz_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_mcz_p_p = kind(mcz_p_p)
  REAL*8 psi (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_psi = kind(psi)
  REAL*8 psix (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_psix = kind(psix)
  REAL*8 psixx (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_psixx = kind(psixx)
  REAL*8 psixy (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_psixy = kind(psixy)
  REAL*8 psixz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_psixz = kind(psixz)
  REAL*8 psiy (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_psiy = kind(psiy)
  REAL*8 psiyy (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_psiyy = kind(psiyy)
  REAL*8 psiyz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_psiyz = kind(psiyz)
  REAL*8 psiz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_psiz = kind(psiz)
  REAL*8 psizz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_psizz = kind(psizz)
  REAL*8 r (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_r = kind(r)
  REAL*8 rhs_alp (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_alp = kind(rhs_alp)
  REAL*8 rhs_axx (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_axx = kind(rhs_axx)
  REAL*8 rhs_axy (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_axy = kind(rhs_axy)
  REAL*8 rhs_axz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_axz = kind(rhs_axz)
  REAL*8 rhs_ayy (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_ayy = kind(rhs_ayy)
  REAL*8 rhs_ayz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_ayz = kind(rhs_ayz)
  REAL*8 rhs_azz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_azz = kind(rhs_azz)
  REAL*8 rhs_betax (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_betax = kind(rhs_betax)
  REAL*8 rhs_betay (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_betay = kind(rhs_betay)
  REAL*8 rhs_betaz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_betaz = kind(rhs_betaz)
  REAL*8 rhs_conf_fac (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_conf_fac = kind(rhs_conf_fac)
  REAL*8 rhs_gammatx (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_gammatx = kind(rhs_gammatx)
  REAL*8 rhs_gammatx_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_gammatx_p = kind(rhs_gammatx_p)
  REAL*8 rhs_gammatx_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_gammatx_p_p = kind(rhs_gammatx_p_p)
  REAL*8 rhs_gammaty (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_gammaty = kind(rhs_gammaty)
  REAL*8 rhs_gammaty_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_gammaty_p = kind(rhs_gammaty_p)
  REAL*8 rhs_gammaty_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_gammaty_p_p = kind(rhs_gammaty_p_p)
  REAL*8 rhs_gammatz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_gammatz = kind(rhs_gammatz)
  REAL*8 rhs_gammatz_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_gammatz_p = kind(rhs_gammatz_p)
  REAL*8 rhs_gammatz_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_gammatz_p_p = kind(rhs_gammatz_p_p)
  REAL*8 rhs_hxx (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_hxx = kind(rhs_hxx)
  REAL*8 rhs_hxy (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_hxy = kind(rhs_hxy)
  REAL*8 rhs_hxz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_hxz = kind(rhs_hxz)
  REAL*8 rhs_hyy (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_hyy = kind(rhs_hyy)
  REAL*8 rhs_hyz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_hyz = kind(rhs_hyz)
  REAL*8 rhs_hzz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_hzz = kind(rhs_hzz)
  REAL*8 rhs_tracek (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_tracek = kind(rhs_tracek)
  INTEGER*4 shift_state
  integer, parameter :: cctki_use_shift_state = kind(shift_state)
  INTEGER*4 stress_energy_state
  integer, parameter :: cctki_use_stress_energy_state = kind(stress_energy_state)
  REAL*8 tracek (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_tracek = kind(tracek)
  REAL*8 tracek_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_tracek_p = kind(tracek_p)
  REAL*8 tracek_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_tracek_p_p = kind(tracek_p_p)
  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)
  INTEGER*4 LeanBSSN_MaxNumConstrainedVars
  integer, parameter :: cctki_use_LeanBSSN_MaxNumConstrainedVars = kind(LeanBSSN_MaxNumConstrainedVars)
  INTEGER*4 LeanBSSN_MaxNumEvolvedVars
  integer, parameter :: cctki_use_LeanBSSN_MaxNumEvolvedVars = kind(LeanBSSN_MaxNumEvolvedVars)
  INTEGER*4 LeanBSSN_MaxNumSandRVars
  integer, parameter :: cctki_use_LeanBSSN_MaxNumSandRVars = kind(LeanBSSN_MaxNumSandRVars)
  COMMON /LeanBSSNMoLrest/LeanBSSN_MaxNumConstrainedVars, LeanBSSN_MaxNumEvolvedVars, LeanBSSN_MaxNumSandRVars
  REAL*8 beta_Alp
  integer, parameter :: cctki_use_beta_Alp = kind(beta_Alp)
  REAL*8 beta_Gamma
  integer, parameter :: cctki_use_beta_Gamma = kind(beta_Gamma)
  REAL*8 beta_f_delta1
  integer, parameter :: cctki_use_beta_f_delta1 = kind(beta_f_delta1)
  REAL*8 beta_f_delta2
  integer, parameter :: cctki_use_beta_f_delta2 = kind(beta_f_delta2)
  REAL*8 chi_gamma
  integer, parameter :: cctki_use_chi_gamma = kind(chi_gamma)
  REAL*8 conf_fac_floor
  integer, parameter :: cctki_use_conf_fac_floor = kind(conf_fac_floor)
  REAL*8 eps_r
  integer, parameter :: cctki_use_eps_r = kind(eps_r)
  REAL*8 eta_beta
  integer, parameter :: cctki_use_eta_beta = kind(eta_beta)
  REAL*8 eta_transition_r
  integer, parameter :: cctki_use_eta_transition_r = kind(eta_transition_r)
  REAL*8 zeta_alpha
  integer, parameter :: cctki_use_zeta_alpha = kind(zeta_alpha)
  REAL*8 zeta_beta
  integer, parameter :: cctki_use_zeta_beta = kind(zeta_beta)
  INTEGER*4 calculate_constraints
  integer, parameter :: cctki_use_calculate_constraints = kind(calculate_constraints)
  INTEGER*4 calculate_constraints_every
  integer, parameter :: cctki_use_calculate_constraints_every = kind(calculate_constraints_every)
  INTEGER*4 compute_rhs_at_initial
  integer, parameter :: cctki_use_compute_rhs_at_initial = kind(compute_rhs_at_initial)
  INTEGER*4 derivs_order
  integer, parameter :: cctki_use_derivs_order = kind(derivs_order)
  INTEGER*4 eta_beta_dynamic
  integer, parameter :: cctki_use_eta_beta_dynamic = kind(eta_beta_dynamic)
  INTEGER*4 eta_transition
  integer, parameter :: cctki_use_eta_transition = kind(eta_transition)
  INTEGER*4 impose_conf_fac_floor_at_initial
  integer, parameter :: cctki_use_impose_conf_fac_floor_at_initial = kind(impose_conf_fac_floor_at_initial)
  INTEGER*4 make_aa_tracefree
  integer, parameter :: cctki_use_make_aa_tracefree = kind(make_aa_tracefree)
  INTEGER*4 moving_eta_transition
  integer, parameter :: cctki_use_moving_eta_transition = kind(moving_eta_transition)
  INTEGER*4 n_aij
  integer, parameter :: cctki_use_n_aij = kind(n_aij)
  INTEGER*4 n_alpha
  integer, parameter :: cctki_use_n_alpha = kind(n_alpha)
  INTEGER*4 n_beta
  integer, parameter :: cctki_use_n_beta = kind(n_beta)
  INTEGER*4 n_conf_fac
  integer, parameter :: cctki_use_n_conf_fac = kind(n_conf_fac)
  INTEGER*4 n_gammat
  integer, parameter :: cctki_use_n_gammat = kind(n_gammat)
  INTEGER*4 n_hij
  integer, parameter :: cctki_use_n_hij = kind(n_hij)
  INTEGER*4 n_trk
  integer, parameter :: cctki_use_n_trk = kind(n_trk)
  INTEGER*4 precollapsed_lapse
  integer, parameter :: cctki_use_precollapsed_lapse = kind(precollapsed_lapse)
  INTEGER*4 rescale_shift_initial
  integer, parameter :: cctki_use_rescale_shift_initial = kind(rescale_shift_initial)
  INTEGER*4 reset_dethh
  integer, parameter :: cctki_use_reset_dethh = kind(reset_dethh)
  INTEGER*4 use_advection_stencils
  integer, parameter :: cctki_use_use_advection_stencils = kind(use_advection_stencils)
  INTEGER*4 z_is_radial
  integer, parameter :: cctki_use_z_is_radial = kind(z_is_radial)
  COMMON /LeanBSSNMoLpriv/beta_Alp, beta_Gamma, beta_f_delta1, beta_f_delta2, chi_gamma, conf_fac_floor, eps_r, eta_beta, eta_trans&
  &ition_r, zeta_alpha, zeta_beta, calculate_constraints, calculate_constraints_every, compute_rhs_at_initial, derivs_order, eta_be&
  &ta_dynamic, eta_transition, impose_conf_fac_floor_at_initial, make_aa_tracefree, moving_eta_transition, n_aij, n_alpha, n_beta, &
  &n_conf_fac, n_gammat, n_hij, n_trk, precollapsed_lapse, rescale_shift_initial, reset_dethh, use_advection_stencils, z_is_radial
  integer*8 CCTKH0
  integer, parameter :: cctki_use_CCTKH0 = kind(CCTKH0)
  integer*8 dtlapse_evolution_method
  integer, parameter :: cctki_use_dtlapse_evolution_method = kind(dtlapse_evolution_method)
  integer*8 dtshift_evolution_method
  integer, parameter :: cctki_use_dtshift_evolution_method = kind(dtshift_evolution_method)
  integer*8 evolution_method
  integer, parameter :: cctki_use_evolution_method = kind(evolution_method)
  integer*8 CCTKH1
  integer, parameter :: cctki_use_CCTKH1 = kind(CCTKH1)
  integer*8 CCTKH2
  integer, parameter :: cctki_use_CCTKH2 = kind(CCTKH2)
  integer*8 CCTKH3
  integer, parameter :: cctki_use_CCTKH3 = kind(CCTKH3)
  integer*8 CCTKH4
  integer, parameter :: cctki_use_CCTKH4 = kind(CCTKH4)
  integer*8 CCTKH5
  integer, parameter :: cctki_use_CCTKH5 = kind(CCTKH5)
  integer*8 lapse_evolution_method
  integer, parameter :: cctki_use_lapse_evolution_method = kind(lapse_evolution_method)
  integer*8 CCTKH6
  integer, parameter :: cctki_use_CCTKH6 = kind(CCTKH6)
  integer*8 CCTKH8
  integer, parameter :: cctki_use_CCTKH8 = kind(CCTKH8)
  integer*8 CCTKH10
  integer, parameter :: cctki_use_CCTKH10 = kind(CCTKH10)
  integer*8 shift_evolution_method
  integer, parameter :: cctki_use_shift_evolution_method = kind(shift_evolution_method)
  integer*8 CCTKH11
  integer, parameter :: cctki_use_CCTKH11 = kind(CCTKH11)
  INTEGER*4 CCTKH7
  integer, parameter :: cctki_use_CCTKH7 = kind(CCTKH7)
  INTEGER*4 CCTKH9
  integer, parameter :: cctki_use_CCTKH9 = kind(CCTKH9)
  INTEGER*4 CCTKH12
  integer, parameter :: cctki_use_CCTKH12 = kind(CCTKH12)
  COMMON /ADMBASErest/CCTKH0, dtlapse_evolution_method, dtshift_evolution_method, evolution_method, CCTKH1, CCTKH2, CCTKH3, CCTKH4,&
  & CCTKH5, lapse_evolution_method, CCTKH6, CCTKH8, CCTKH10, shift_evolution_method, CCTKH11, CCTKH7, CCTKH9, CCTKH12
  INTEGER*4 CCTKH13
  integer, parameter :: cctki_use_CCTKH13 = kind(CCTKH13)
  INTEGER*4 CCTKH14
  integer, parameter :: cctki_use_CCTKH14 = kind(CCTKH14)
  INTEGER*4 CCTKH15
  integer, parameter :: cctki_use_CCTKH15 = kind(CCTKH15)
  INTEGER*4 CCTKH16
  integer, parameter :: cctki_use_CCTKH16 = kind(CCTKH16)
  INTEGER*4 MoL_Num_Constrained_Vars
  integer, parameter :: cctki_use_MoL_Num_Constrained_Vars = kind(MoL_Num_Constrained_Vars)
  INTEGER*4 MoL_Num_Evolved_Vars
  integer, parameter :: cctki_use_MoL_Num_Evolved_Vars = kind(MoL_Num_Evolved_Vars)
  INTEGER*4 CCTKH17
  integer, parameter :: cctki_use_CCTKH17 = kind(CCTKH17)
  INTEGER*4 MoL_Num_SaveAndRestore_Vars
  integer, parameter :: cctki_use_MoL_Num_SaveAndRestore_Vars = kind(MoL_Num_SaveAndRestore_Vars)
  INTEGER*4 CCTKH18
  integer, parameter :: cctki_use_CCTKH18 = kind(CCTKH18)
  COMMON /METHODOFLINESrest/CCTKH13, CCTKH14, CCTKH15, CCTKH16, MoL_Num_Constrained_Vars, MoL_Num_Evolved_Vars, CCTKH17, MoL_Num_Sa&
  &veAndRestore_Vars, CCTKH18
  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
  interface
  INTEGER*4 function ExtrapolateGammas (cctkGH, var)
  implicit none
  integer*8 cctkGH
  REAL*8 var(*)
  end function ExtrapolateGammas
  end interface
  interface
  INTEGER*4 function GetRefinementLevel (cctkGH)
  implicit none
  integer*8 cctkGH
  end function GetRefinementLevel
  end interface
  interface
  INTEGER*4 function MoLRegisterConstrainedGroup (ConstrainedIndex)
  implicit none
  INTEGER*4 ConstrainedIndex
  end function MoLRegisterConstrainedGroup
  end interface
  interface
  INTEGER*4 function MoLRegisterEvolved (EvolvedIndex, RHSIndex)
  implicit none
  INTEGER*4 EvolvedIndex
  INTEGER*4 RHSIndex
  end function MoLRegisterEvolved
  end interface
  interface
  INTEGER*4 function MoLRegisterEvolvedGroup (EvolvedIndex, RHSIndex)
  implicit none
  INTEGER*4 EvolvedIndex
  INTEGER*4 RHSIndex
  end function MoLRegisterEvolvedGroup
  end interface
  interface
  INTEGER*4 function MoLRegisterSaveAndRestoreGroup (SandRIndex)
  implicit none
  INTEGER*4 SandRIndex
  end function MoLRegisterSaveAndRestoreGroup
  end interface
  interface
  INTEGER*4 function MultiPatch_GetMap (cctkGH)
  implicit none
  integer*8 cctkGH
  end function MultiPatch_GetMap
  end interface
  interface
  INTEGER*4 function NewRad_Apply (cctkGH, var, rhs, var0, v0, radpower)
  implicit none
  integer*8 cctkGH
  REAL*8 var(*)
  REAL*8 rhs(*)
  REAL*8 var0
  REAL*8 v0
  INTEGER*4 radpower
  end function NewRad_Apply
  end interface

  
  REAL*8                alph, beta(3)
  REAL*8                ww, hh(3,3), hu(3,3), trk, aa(3,3), gammat(3),   &
                           dethh, Tab(4,4)

  
  REAL*8                d1_hh11(3), d1_hh12(3), d1_hh13(3), d1_hh22(3), d1_hh23(3), d1_hh33(3)
  REAL*8                d1_aa11(3), d1_aa12(3), d1_aa13(3), d1_aa22(3), d1_aa23(3), d1_aa33(3)
  REAL*8                d1_gammat1(3), d1_gammat2(3), d1_gammat3(3)
  REAL*8                d1_ww(3), d1_hh(3,3,3), d1_trk(3), d1_aa(3,3,3),&
                           d1_gammat(3,3)

  
  REAL*8                d2_hh11(3,3), d2_hh12(3,3), d2_hh13(3,3), d2_hh22(3,3), d2_hh23(3,3), d2_hh33(3,3)
  REAL*8                d2_ww(3,3), d2_hh(3,3,3,3)


  
  REAL*8                cd2_ww(3,3), cd1_aa(3,3,3)

  
  REAL*8                ham, mom(3)

  
  REAL*8                cf1(3,3,3), cf2(3,3,3), c_ri(3,3)
  REAL*8                c_ri_ww(3,3), c_ri_hh(3,3), ri_1(3,3), ri_2(3,3),&
                           ri_3(3,3), sq_aa, a2(3,3), trr
  REAL*8                tr_cd2_ww, tr_dww_dww

  
  REAL*8                srcE, srcS, srcjdi(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, parameter ::  one  = 1
  REAL*8, parameter ::  pi   = acos(-one)
  REAL*8, parameter ::  pi8  = 8*pi
  REAL*8, parameter ::  pi16 = 16*pi
  INTEGER*4                 i, j, k
  INTEGER*4                 a, b, c, l, m, n, p, q

  
  REAL*8                jac(3,3), hes(3,3,3)

  integer                  istat
  logical                  use_jacobian
  REAL*8, dimension(cctk_ash(1),cctk_ash(2),cctk_ash(3)) :: lJ11, lJ12, lJ13
  REAL*8, dimension(cctk_ash(1),cctk_ash(2),cctk_ash(3)) :: lJ21, lJ22, lJ23
  REAL*8, dimension(cctk_ash(1),cctk_ash(2),cctk_ash(3)) :: lJ31, lJ32, lJ33
  REAL*8, dimension(cctk_ash(1),cctk_ash(2),cctk_ash(3)) :: ldJ111, ldJ112, ldJ113, ldJ122, ldJ123, ldJ133
  REAL*8, dimension(cctk_ash(1),cctk_ash(2),cctk_ash(3)) :: ldJ211, ldJ212, ldJ213, ldJ222, ldJ223, ldJ233
  REAL*8, dimension(cctk_ash(1),cctk_ash(2),cctk_ash(3)) :: ldJ311, ldJ312, ldJ313, ldJ322, ldJ323, ldJ333

  integer*8             lJ11_ptr, lJ12_ptr, lJ13_ptr
  integer*8             lJ21_ptr, lJ22_ptr, lJ23_ptr
  integer*8             lJ31_ptr, lJ32_ptr, lJ33_ptr
  integer*8             ldJ111_ptr, ldJ112_ptr, ldJ113_ptr, ldJ122_ptr, ldJ123_ptr, ldJ133_ptr
  integer*8             ldJ211_ptr, ldJ212_ptr, ldJ213_ptr, ldJ222_ptr, ldJ223_ptr, ldJ233_ptr
  integer*8             ldJ311_ptr, ldJ312_ptr, ldJ313_ptr, ldJ322_ptr, ldJ323_ptr, ldJ333_ptr

  pointer (lJ11_ptr, lJ11), (lJ12_ptr, lJ12), (lJ13_ptr, lJ13)
  pointer (lJ21_ptr, lJ21), (lJ22_ptr, lJ22), (lJ23_ptr, lJ23)
  pointer (lJ31_ptr, lJ31), (lJ32_ptr, lJ32), (lJ33_ptr, lJ33)

  pointer (ldJ111_ptr, ldJ111), (ldJ112_ptr, ldJ112), (ldJ113_ptr, ldJ113), (ldJ122_ptr, ldJ122), (ldJ123_ptr, ldJ123), (ldJ133_ptr&
  &, ldJ133)
  pointer (ldJ211_ptr, ldJ211), (ldJ212_ptr, ldJ212), (ldJ213_ptr, ldJ213), (ldJ222_ptr, ldJ222), (ldJ223_ptr, ldJ223), (ldJ233_ptr&
  &, ldJ233)
  pointer (ldJ311_ptr, ldJ311), (ldJ312_ptr, ldJ312), (ldJ313_ptr, ldJ313), (ldJ322_ptr, ldJ322), (ldJ323_ptr, ldJ323), (ldJ333_ptr&
  &, ldJ333)

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

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

  
  call CCTK_IsFunctionAliased(istat, "MultiPatch_GetDomainSpecification")
  if (istat == 0) then
     use_jacobian = .false.
  else
     use_jacobian = .true.
  end if

  if (use_jacobian) then
     call CCTK_VarDataPtr(lJ11_ptr, cctkGH, 0, "Coordinates::J11")
     call CCTK_VarDataPtr(lJ12_ptr, cctkGH, 0, "Coordinates::J12")
     call CCTK_VarDataPtr(lJ13_ptr, cctkGH, 0, "Coordinates::J13")
     call CCTK_VarDataPtr(lJ21_ptr, cctkGH, 0, "Coordinates::J21")
     call CCTK_VarDataPtr(lJ22_ptr, cctkGH, 0, "Coordinates::J22")
     call CCTK_VarDataPtr(lJ23_ptr, cctkGH, 0, "Coordinates::J23")
     call CCTK_VarDataPtr(lJ31_ptr, cctkGH, 0, "Coordinates::J31")
     call CCTK_VarDataPtr(lJ32_ptr, cctkGH, 0, "Coordinates::J32")
     call CCTK_VarDataPtr(lJ33_ptr, cctkGH, 0, "Coordinates::J33")

     call CCTK_VarDataPtr(ldJ111_ptr, cctkGH, 0, "Coordinates::dJ111")
     call CCTK_VarDataPtr(ldJ112_ptr, cctkGH, 0, "Coordinates::dJ112")
     call CCTK_VarDataPtr(ldJ113_ptr, cctkGH, 0, "Coordinates::dJ113")
     call CCTK_VarDataPtr(ldJ122_ptr, cctkGH, 0, "Coordinates::dJ122")
     call CCTK_VarDataPtr(ldJ123_ptr, cctkGH, 0, "Coordinates::dJ123")
     call CCTK_VarDataPtr(ldJ133_ptr, cctkGH, 0, "Coordinates::dJ133")

     call CCTK_VarDataPtr(ldJ211_ptr, cctkGH, 0, "Coordinates::dJ211")
     call CCTK_VarDataPtr(ldJ212_ptr, cctkGH, 0, "Coordinates::dJ212")
     call CCTK_VarDataPtr(ldJ213_ptr, cctkGH, 0, "Coordinates::dJ213")
     call CCTK_VarDataPtr(ldJ222_ptr, cctkGH, 0, "Coordinates::dJ222")
     call CCTK_VarDataPtr(ldJ223_ptr, cctkGH, 0, "Coordinates::dJ223")
     call CCTK_VarDataPtr(ldJ233_ptr, cctkGH, 0, "Coordinates::dJ233")

     call CCTK_VarDataPtr(ldJ311_ptr, cctkGH, 0, "Coordinates::dJ311")
     call CCTK_VarDataPtr(ldJ312_ptr, cctkGH, 0, "Coordinates::dJ312")
     call CCTK_VarDataPtr(ldJ313_ptr, cctkGH, 0, "Coordinates::dJ313")
     call CCTK_VarDataPtr(ldJ322_ptr, cctkGH, 0, "Coordinates::dJ322")
     call CCTK_VarDataPtr(ldJ323_ptr, cctkGH, 0, "Coordinates::dJ323")
     call CCTK_VarDataPtr(ldJ333_ptr, cctkGH, 0, "Coordinates::dJ333")
  end if

  dx12 = 12*(cctk_delta_space(1)/cctk_levfac(1))
  dy12 = 12*(cctk_delta_space(2)/cctk_levfac(2))
  dz12 = 12*(cctk_delta_space(3)/cctk_levfac(3))

  dxsq12 = 12*(cctk_delta_space(1)/cctk_levfac(1))*(cctk_delta_space(1)/cctk_levfac(1))
  dysq12 = 12*(cctk_delta_space(2)/cctk_levfac(2))*(cctk_delta_space(2)/cctk_levfac(2))
  dzsq12 = 12*(cctk_delta_space(3)/cctk_levfac(3))*(cctk_delta_space(3)/cctk_levfac(3))

  dxdy144 = 144*(cctk_delta_space(1)/cctk_levfac(1))*(cctk_delta_space(2)/cctk_levfac(2))
  dxdz144 = 144*(cctk_delta_space(1)/cctk_levfac(1))*(cctk_delta_space(3)/cctk_levfac(3))
  dydz144 = 144*(cctk_delta_space(2)/cctk_levfac(2))*(cctk_delta_space(3)/cctk_levfac(3))


  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)))

  
  hc    = 0
  mcx   = 0
  mcy   = 0
  mcz   = 0

  !$OMP PARALLEL DO COLLAPSE(3) &
  !$OMP PRIVATE( alph, beta, &
  !$OMP ww, hh, hu, trk, aa, gammat, dethh, Tab, &
  !$OMP d1_hh11, d1_hh12, d1_hh13, d1_hh22, d1_hh23, d1_hh33, &
  !$OMP d1_aa11, d1_aa12, d1_aa13, d1_aa22, d1_aa23, d1_aa33, &
  !$OMP d1_gammat1, d1_gammat2, d1_gammat3, &
  !$OMP d1_ww, d1_hh, d1_trk, d1_aa, d1_gammat, &
  !$OMP d2_ww, d2_hh, cd2_ww, cd1_aa, ham, mom, &
  !$OMP cf1, cf2, c_ri, c_ri_ww, c_ri_hh, ri_1, &
  !$OMP ri_2, ri_3, sq_aa, a2, trr, &
  !$OMP tr_cd2_ww, tr_dww_dww, srcE, srcS, srcjdi, &
  !$OMP i, j, k, a, b, c, l, m, n, p, q, jac, hes)
  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)

    
    ww        = conf_fac(i,j,k)

    hh(1,1)   = hxx(i,j,k)
    hh(1,2)   = hxy(i,j,k)
    hh(1,3)   = hxz(i,j,k)
    hh(2,2)   = hyy(i,j,k)
    hh(2,3)   = hyz(i,j,k)
    hh(3,3)   = hzz(i,j,k)
    hh(2,1)   = hh(1,2)
    hh(3,1)   = hh(1,3)
    hh(3,2)   = hh(2,3)

    trk       = tracek(i,j,k)

    aa(1,1)   = axx(i,j,k)
    aa(1,2)   = axy(i,j,k)
    aa(1,3)   = axz(i,j,k)
    aa(2,2)   = ayy(i,j,k)
    aa(2,3)   = ayz(i,j,k)
    aa(3,3)   = azz(i,j,k)
    aa(2,1)   = aa(1,2)
    aa(3,1)   = aa(1,3)
    aa(3,2)   = aa(2,3)

    gammat(1) = gammatx(i,j,k)
    gammat(2) = gammaty(i,j,k)
    gammat(3) = gammatz(i,j,k)

    alph      = alp(i,j,k)

    beta(1)   = betax(i,j,k)
    beta(2)   = betay(i,j,k)
    beta(3)   = betaz(i,j,k)

    if (use_jacobian) then
       jac(1,1) = lJ11(i,j,k)
       jac(1,2) = lJ12(i,j,k)
       jac(1,3) = lJ13(i,j,k)
       jac(2,1) = lJ21(i,j,k)
       jac(2,2) = lJ22(i,j,k)
       jac(2,3) = lJ23(i,j,k)
       jac(3,1) = lJ31(i,j,k)
       jac(3,2) = lJ32(i,j,k)
       jac(3,3) = lJ33(i,j,k)

       hes(1,1,1) = ldJ111(i,j,k)
       hes(1,1,2) = ldJ112(i,j,k)
       hes(1,1,3) = ldJ113(i,j,k)
       hes(1,2,1) = ldJ112(i,j,k)
       hes(1,2,2) = ldJ122(i,j,k)
       hes(1,2,3) = ldJ123(i,j,k)
       hes(1,3,1) = ldJ113(i,j,k)
       hes(1,3,2) = ldJ123(i,j,k)
       hes(1,3,3) = ldJ133(i,j,k)

       hes(2,1,1) = ldJ211(i,j,k)
       hes(2,1,2) = ldJ212(i,j,k)
       hes(2,1,3) = ldJ213(i,j,k)
       hes(2,2,1) = ldJ212(i,j,k)
       hes(2,2,2) = ldJ222(i,j,k)
       hes(2,2,3) = ldJ223(i,j,k)
       hes(2,3,1) = ldJ213(i,j,k)
       hes(2,3,2) = ldJ223(i,j,k)
       hes(2,3,3) = ldJ233(i,j,k)

       hes(3,1,1) = ldJ311(i,j,k)
       hes(3,1,2) = ldJ312(i,j,k)
       hes(3,1,3) = ldJ313(i,j,k)
       hes(3,2,1) = ldJ312(i,j,k)
       hes(3,2,2) = ldJ322(i,j,k)
       hes(3,2,3) = ldJ323(i,j,k)
       hes(3,3,1) = ldJ313(i,j,k)
       hes(3,3,2) = ldJ323(i,j,k)
       hes(3,3,3) = ldJ333(i,j,k)
    else
       jac      = 0.0
       jac(1,1) = 1.0
       jac(2,2) = 1.0
       jac(3,3) = 1.0
       hes      = 0.0
    end if

    
    
    dethh =       hh(1,1) * hh(2,2) * hh(3,3)                              &
            + 2 * hh(1,2) * hh(1,3) * hh(2,3)                              &
            -     hh(1,1) * hh(2,3) ** 2                                   &
            -     hh(2,2) * hh(1,3) ** 2                                   &
            -     hh(3,3) * hh(1,2) ** 2
    hu(1,1) = (hh(2,2) * hh(3,3) - hh(2,3) ** 2     ) / dethh
    hu(2,2) = (hh(1,1) * hh(3,3) - hh(1,3) ** 2     ) / dethh
    hu(3,3) = (hh(1,1) * hh(2,2) - hh(1,2) ** 2     ) / dethh
    hu(1,2) = (hh(1,3) * hh(2,3) - hh(1,2) * hh(3,3)) / dethh
    hu(1,3) = (hh(1,2) * hh(2,3) - hh(1,3) * hh(2,2)) / dethh
    hu(2,3) = (hh(1,3) * hh(1,2) - hh(2,3) * hh(1,1)) / dethh
    hu(2,1) = hu(1,2)
    hu(3,1) = hu(1,3)
    hu(3,2) = hu(2,3)
    


    if (derivs_order == 6) then

      
      
      d1_ww(1) = (  conf_fac(i+3,j,k) - 9*conf_fac(i+2,j,k) + 45*conf_fac(i+1,j,k)          &
                  - conf_fac(i-3,j,k) + 9*conf_fac(i-2,j,k) - 45*conf_fac(i-1,j,k) ) * odx60

      d1_ww(2) = (  conf_fac(i,j+3,k) - 9*conf_fac(i,j+2,k) + 45*conf_fac(i,j+1,k)          &
                  - conf_fac(i,j-3,k) + 9*conf_fac(i,j-2,k) - 45*conf_fac(i,j-1,k) ) * ody60

      d1_ww(3) = (  conf_fac(i,j,k+3) - 9*conf_fac(i,j,k+2) + 45*conf_fac(i,j,k+1)          &
                  - conf_fac(i,j,k-3) + 9*conf_fac(i,j,k-2) - 45*conf_fac(i,j,k-1) ) * odz60

      
      d1_hh11(1) = (  hxx(i+3,j,k) - 9*hxx(i+2,j,k) + 45*hxx(i+1,j,k)      &
                    - hxx(i-3,j,k) + 9*hxx(i-2,j,k) - 45*hxx(i-1,j,k) ) * odx60
      d1_hh12(1) = (  hxy(i+3,j,k) - 9*hxy(i+2,j,k) + 45*hxy(i+1,j,k)      &
                    - hxy(i-3,j,k) + 9*hxy(i-2,j,k) - 45*hxy(i-1,j,k) ) * odx60
      d1_hh13(1) = (  hxz(i+3,j,k) - 9*hxz(i+2,j,k) + 45*hxz(i+1,j,k)      &
                    - hxz(i-3,j,k) + 9*hxz(i-2,j,k) - 45*hxz(i-1,j,k) ) * odx60
      d1_hh22(1) = (  hyy(i+3,j,k) - 9*hyy(i+2,j,k) + 45*hyy(i+1,j,k)      &
                    - hyy(i-3,j,k) + 9*hyy(i-2,j,k) - 45*hyy(i-1,j,k) ) * odx60
      d1_hh23(1) = (  hyz(i+3,j,k) - 9*hyz(i+2,j,k) + 45*hyz(i+1,j,k)      &
                    - hyz(i-3,j,k) + 9*hyz(i-2,j,k) - 45*hyz(i-1,j,k) ) * odx60
      d1_hh33(1) = (  hzz(i+3,j,k) - 9*hzz(i+2,j,k) + 45*hzz(i+1,j,k)      &
                    - hzz(i-3,j,k) + 9*hzz(i-2,j,k) - 45*hzz(i-1,j,k) ) * odx60

      d1_hh11(2) = (  hxx(i,j+3,k) - 9*hxx(i,j+2,k) + 45*hxx(i,j+1,k)      &
                    - hxx(i,j-3,k) + 9*hxx(i,j-2,k) - 45*hxx(i,j-1,k) ) * ody60
      d1_hh12(2) = (  hxy(i,j+3,k) - 9*hxy(i,j+2,k) + 45*hxy(i,j+1,k)      &
                    - hxy(i,j-3,k) + 9*hxy(i,j-2,k) - 45*hxy(i,j-1,k) ) * ody60
      d1_hh13(2) = (  hxz(i,j+3,k) - 9*hxz(i,j+2,k) + 45*hxz(i,j+1,k)      &
                    - hxz(i,j-3,k) + 9*hxz(i,j-2,k) - 45*hxz(i,j-1,k) ) * ody60
      d1_hh22(2) = (  hyy(i,j+3,k) - 9*hyy(i,j+2,k) + 45*hyy(i,j+1,k)      &
                    - hyy(i,j-3,k) + 9*hyy(i,j-2,k) - 45*hyy(i,j-1,k) ) * ody60
      d1_hh23(2) = (  hyz(i,j+3,k) - 9*hyz(i,j+2,k) + 45*hyz(i,j+1,k)      &
                    - hyz(i,j-3,k) + 9*hyz(i,j-2,k) - 45*hyz(i,j-1,k) ) * ody60
      d1_hh33(2) = (  hzz(i,j+3,k) - 9*hzz(i,j+2,k) + 45*hzz(i,j+1,k)      &
                    - hzz(i,j-3,k) + 9*hzz(i,j-2,k) - 45*hzz(i,j-1,k) ) * ody60

      d1_hh11(3) = (  hxx(i,j,k+3) - 9*hxx(i,j,k+2) + 45*hxx(i,j,k+1)      &
                    - hxx(i,j,k-3) + 9*hxx(i,j,k-2) - 45*hxx(i,j,k-1) ) * odz60
      d1_hh12(3) = (  hxy(i,j,k+3) - 9*hxy(i,j,k+2) + 45*hxy(i,j,k+1)      &
                    - hxy(i,j,k-3) + 9*hxy(i,j,k-2) - 45*hxy(i,j,k-1) ) * odz60
      d1_hh13(3) = (  hxz(i,j,k+3) - 9*hxz(i,j,k+2) + 45*hxz(i,j,k+1)      &
                    - hxz(i,j,k-3) + 9*hxz(i,j,k-2) - 45*hxz(i,j,k-1) ) * odz60
      d1_hh22(3) = (  hyy(i,j,k+3) - 9*hyy(i,j,k+2) + 45*hyy(i,j,k+1)      &
                    - hyy(i,j,k-3) + 9*hyy(i,j,k-2) - 45*hyy(i,j,k-1) ) * odz60
      d1_hh23(3) = (  hyz(i,j,k+3) - 9*hyz(i,j,k+2) + 45*hyz(i,j,k+1)      &
                    - hyz(i,j,k-3) + 9*hyz(i,j,k-2) - 45*hyz(i,j,k-1) ) * odz60
      d1_hh33(3) = (  hzz(i,j,k+3) - 9*hzz(i,j,k+2) + 45*hzz(i,j,k+1)      &
                    - hzz(i,j,k-3) + 9*hzz(i,j,k-2) - 45*hzz(i,j,k-1) ) * odz60

      
      d1_trk(1) = (  tracek(i+3,j,k) - 9*tracek(i+2,j,k) + 45*tracek(i+1,j,k)          &
                   - tracek(i-3,j,k) + 9*tracek(i-2,j,k) - 45*tracek(i-1,j,k) ) * odx60

      d1_trk(2) = (  tracek(i,j+3,k) - 9*tracek(i,j+2,k) + 45*tracek(i,j+1,k)          &
                   - tracek(i,j-3,k) + 9*tracek(i,j-2,k) - 45*tracek(i,j-1,k) ) * ody60

      d1_trk(3) = (  tracek(i,j,k+3) - 9*tracek(i,j,k+2) + 45*tracek(i,j,k+1)          &
                   - tracek(i,j,k-3) + 9*tracek(i,j,k-2) - 45*tracek(i,j,k-1) ) * odz60


      
      d1_aa11(1) = (  axx(i+3,j,k) - 9*axx(i+2,j,k) + 45*axx(i+1,j,k)      &
                    - axx(i-3,j,k) + 9*axx(i-2,j,k) - 45*axx(i-1,j,k) ) * odx60
      d1_aa12(1) = (  axy(i+3,j,k) - 9*axy(i+2,j,k) + 45*axy(i+1,j,k)      &
                    - axy(i-3,j,k) + 9*axy(i-2,j,k) - 45*axy(i-1,j,k) ) * odx60
      d1_aa13(1) = (  axz(i+3,j,k) - 9*axz(i+2,j,k) + 45*axz(i+1,j,k)      &
                    - axz(i-3,j,k) + 9*axz(i-2,j,k) - 45*axz(i-1,j,k) ) * odx60
      d1_aa22(1) = (  ayy(i+3,j,k) - 9*ayy(i+2,j,k) + 45*ayy(i+1,j,k)      &
                    - ayy(i-3,j,k) + 9*ayy(i-2,j,k) - 45*ayy(i-1,j,k) ) * odx60
      d1_aa23(1) = (  ayz(i+3,j,k) - 9*ayz(i+2,j,k) + 45*ayz(i+1,j,k)      &
                    - ayz(i-3,j,k) + 9*ayz(i-2,j,k) - 45*ayz(i-1,j,k) ) * odx60
      d1_aa33(1) = (  azz(i+3,j,k) - 9*azz(i+2,j,k) + 45*azz(i+1,j,k)      &
                    - azz(i-3,j,k) + 9*azz(i-2,j,k) - 45*azz(i-1,j,k) ) * odx60

      d1_aa11(2) = (  axx(i,j+3,k) - 9*axx(i,j+2,k) + 45*axx(i,j+1,k)      &
                    - axx(i,j-3,k) + 9*axx(i,j-2,k) - 45*axx(i,j-1,k) ) * ody60
      d1_aa12(2) = (  axy(i,j+3,k) - 9*axy(i,j+2,k) + 45*axy(i,j+1,k)      &
                    - axy(i,j-3,k) + 9*axy(i,j-2,k) - 45*axy(i,j-1,k) ) * ody60
      d1_aa13(2) = (  axz(i,j+3,k) - 9*axz(i,j+2,k) + 45*axz(i,j+1,k)      &
                    - axz(i,j-3,k) + 9*axz(i,j-2,k) - 45*axz(i,j-1,k) ) * ody60
      d1_aa22(2) = (  ayy(i,j+3,k) - 9*ayy(i,j+2,k) + 45*ayy(i,j+1,k)      &
                    - ayy(i,j-3,k) + 9*ayy(i,j-2,k) - 45*ayy(i,j-1,k) ) * ody60
      d1_aa23(2) = (  ayz(i,j+3,k) - 9*ayz(i,j+2,k) + 45*ayz(i,j+1,k)      &
                    - ayz(i,j-3,k) + 9*ayz(i,j-2,k) - 45*ayz(i,j-1,k) ) * ody60
      d1_aa33(2) = (  azz(i,j+3,k) - 9*azz(i,j+2,k) + 45*azz(i,j+1,k)      &
                    - azz(i,j-3,k) + 9*azz(i,j-2,k) - 45*azz(i,j-1,k) ) * ody60

      d1_aa11(3) = (  axx(i,j,k+3) - 9*axx(i,j,k+2) + 45*axx(i,j,k+1)      &
                    - axx(i,j,k-3) + 9*axx(i,j,k-2) - 45*axx(i,j,k-1) ) * odz60
      d1_aa12(3) = (  axy(i,j,k+3) - 9*axy(i,j,k+2) + 45*axy(i,j,k+1)      &
                    - axy(i,j,k-3) + 9*axy(i,j,k-2) - 45*axy(i,j,k-1) ) * odz60
      d1_aa13(3) = (  axz(i,j,k+3) - 9*axz(i,j,k+2) + 45*axz(i,j,k+1)      &
                    - axz(i,j,k-3) + 9*axz(i,j,k-2) - 45*axz(i,j,k-1) ) * odz60
      d1_aa22(3) = (  ayy(i,j,k+3) - 9*ayy(i,j,k+2) + 45*ayy(i,j,k+1)      &
                    - ayy(i,j,k-3) + 9*ayy(i,j,k-2) - 45*ayy(i,j,k-1) ) * odz60
      d1_aa23(3) = (  ayz(i,j,k+3) - 9*ayz(i,j,k+2) + 45*ayz(i,j,k+1)      &
                    - ayz(i,j,k-3) + 9*ayz(i,j,k-2) - 45*ayz(i,j,k-1) ) * odz60
      d1_aa33(3) = (  azz(i,j,k+3) - 9*azz(i,j,k+2) + 45*azz(i,j,k+1)      &
                    - azz(i,j,k-3) + 9*azz(i,j,k-2) - 45*azz(i,j,k-1) ) * odz60

      
      d1_gammat1(1) = (  gammatx(i+3,j,k) - 9*gammatx(i+2,j,k) + 45*gammatx(i+1,j,k) &
                       - gammatx(i-3,j,k) + 9*gammatx(i-2,j,k) - 45*gammatx(i-1,j,k) ) * odx60
      d1_gammat2(1) = (  gammaty(i+3,j,k) - 9*gammaty(i+2,j,k) + 45*gammaty(i+1,j,k) &
                       - gammaty(i-3,j,k) + 9*gammaty(i-2,j,k) - 45*gammaty(i-1,j,k) ) * odx60
      d1_gammat3(1) = (  gammatz(i+3,j,k) - 9*gammatz(i+2,j,k) + 45*gammatz(i+1,j,k) &
                       - gammatz(i-3,j,k) + 9*gammatz(i-2,j,k) - 45*gammatz(i-1,j,k) ) * odx60

      d1_gammat1(2) = (  gammatx(i,j+3,k) - 9*gammatx(i,j+2,k) + 45*gammatx(i,j+1,k) &
                       - gammatx(i,j-3,k) + 9*gammatx(i,j-2,k) - 45*gammatx(i,j-1,k) ) * ody60
      d1_gammat2(2) = (  gammaty(i,j+3,k) - 9*gammaty(i,j+2,k) + 45*gammaty(i,j+1,k) &
                       - gammaty(i,j-3,k) + 9*gammaty(i,j-2,k) - 45*gammaty(i,j-1,k) ) * ody60
      d1_gammat3(2) = (  gammatz(i,j+3,k) - 9*gammatz(i,j+2,k) + 45*gammatz(i,j+1,k) &
                       - gammatz(i,j-3,k) + 9*gammatz(i,j-2,k) - 45*gammatz(i,j-1,k) ) * ody60

      d1_gammat1(3) = (  gammatx(i,j,k+3) - 9*gammatx(i,j,k+2) + 45*gammatx(i,j,k+1) &
                       - gammatx(i,j,k-3) + 9*gammatx(i,j,k-2) - 45*gammatx(i,j,k-1) ) * odz60
      d1_gammat2(3) = (  gammaty(i,j,k+3) - 9*gammaty(i,j,k+2) + 45*gammaty(i,j,k+1) &
                       - gammaty(i,j,k-3) + 9*gammaty(i,j,k-2) - 45*gammaty(i,j,k-1) ) * odz60
      d1_gammat3(3) = (  gammatz(i,j,k+3) - 9*gammatz(i,j,k+2) + 45*gammatz(i,j,k+1) &
                       - gammatz(i,j,k-3) + 9*gammatz(i,j,k-2) - 45*gammatz(i,j,k-1) ) * odz60

      


      

      
      d2_ww(1,1) = (  2*conf_fac(i+3,j,k) - 27*conf_fac(i+2,j,k) + 270*conf_fac(i+1,j,k) - 490*conf_fac(i,j,k)&
                    + 2*conf_fac(i-3,j,k) - 27*conf_fac(i-2,j,k) + 270*conf_fac(i-1,j,k) ) * odxsq180

      d2_ww(2,2) = (  2*conf_fac(i,j+3,k) - 27*conf_fac(i,j+2,k) + 270*conf_fac(i,j+1,k) - 490*conf_fac(i,j,k)&
                    + 2*conf_fac(i,j-3,k) - 27*conf_fac(i,j-2,k) + 270*conf_fac(i,j-1,k) ) * odysq180

      d2_ww(3,3) = (  2*conf_fac(i,j,k+3) - 27*conf_fac(i,j,k+2) + 270*conf_fac(i,j,k+1) - 490*conf_fac(i,j,k)&
                    + 2*conf_fac(i,j,k-3) - 27*conf_fac(i,j,k-2) + 270*conf_fac(i,j,k-1) ) * odzsq180

      d2_ww(1,2) = (    -conf_fac(i-3,j+3,k) +   9*conf_fac(i-2,j+3,k) -   45*conf_fac(i-1,j+3,k) +   45*conf_fac(i+1,j+3,k) -   9*&
  &conf_fac(i+2,j+3,k) +    conf_fac(i+3,j+3,k) &
                    +  9*conf_fac(i-3,j+2,k) -  81*conf_fac(i-2,j+2,k) +  405*conf_fac(i-1,j+2,k) -  405*conf_fac(i+1,j+2,k) +  81*&
  &conf_fac(i+2,j+2,k) -  9*conf_fac(i+3,j+2,k) &
                    - 45*conf_fac(i-3,j+1,k) + 405*conf_fac(i-2,j+1,k) - 2025*conf_fac(i-1,j+1,k) + 2025*conf_fac(i+1,j+1,k) - 405*&
  &conf_fac(i+2,j+1,k) + 45*conf_fac(i+3,j+1,k) &
                    + 45*conf_fac(i-3,j-1,k) - 405*conf_fac(i-2,j-1,k) + 2025*conf_fac(i-1,j-1,k) - 2025*conf_fac(i+1,j-1,k) + 405*&
  &conf_fac(i+2,j-1,k) - 45*conf_fac(i+3,j-1,k) &
                    -  9*conf_fac(i-3,j-2,k) +  81*conf_fac(i-2,j-2,k) -  405*conf_fac(i-1,j-2,k) +  405*conf_fac(i+1,j-2,k) -  81*&
  &conf_fac(i+2,j-2,k) +  9*conf_fac(i+3,j-2,k) &
                    +    conf_fac(i-3,j-3,k) -   9*conf_fac(i-2,j-3,k) +   45*conf_fac(i-1,j-3,k) -   45*conf_fac(i+1,j-3,k) +   9*&
  &conf_fac(i+2,j-3,k) -    conf_fac(i+3,j-3,k) ) * odxdy3600

      d2_ww(1,3) = (    -conf_fac(i-3,j,k+3) +   9*conf_fac(i-2,j,k+3) -   45*conf_fac(i-1,j,k+3) +   45*conf_fac(i+1,j,k+3) -   9*&
  &conf_fac(i+2,j,k+3) +    conf_fac(i+3,j,k+3) &
                    +  9*conf_fac(i-3,j,k+2) -  81*conf_fac(i-2,j,k+2) +  405*conf_fac(i-1,j,k+2) -  405*conf_fac(i+1,j,k+2) +  81*&
  &conf_fac(i+2,j,k+2) -  9*conf_fac(i+3,j,k+2) &
                    - 45*conf_fac(i-3,j,k+1) + 405*conf_fac(i-2,j,k+1) - 2025*conf_fac(i-1,j,k+1) + 2025*conf_fac(i+1,j,k+1) - 405*&
  &conf_fac(i+2,j,k+1) + 45*conf_fac(i+3,j,k+1) &
                    + 45*conf_fac(i-3,j,k-1) - 405*conf_fac(i-2,j,k-1) + 2025*conf_fac(i-1,j,k-1) - 2025*conf_fac(i+1,j,k-1) + 405*&
  &conf_fac(i+2,j,k-1) - 45*conf_fac(i+3,j,k-1) &
                    -  9*conf_fac(i-3,j,k-2) +  81*conf_fac(i-2,j,k-2) -  405*conf_fac(i-1,j,k-2) +  405*conf_fac(i+1,j,k-2) -  81*&
  &conf_fac(i+2,j,k-2) +  9*conf_fac(i+3,j,k-2) &
                    +    conf_fac(i-3,j,k-3) -   9*conf_fac(i-2,j,k-3) +   45*conf_fac(i-1,j,k-3) -   45*conf_fac(i+1,j,k-3) +   9*&
  &conf_fac(i+2,j,k-3) -    conf_fac(i+3,j,k-3) ) * odxdz3600

      d2_ww(2,3) = (    -conf_fac(i,j-3,k+3) +   9*conf_fac(i,j-2,k+3) -   45*conf_fac(i,j-1,k+3) +   45*conf_fac(i,j+1,k+3) -   9*&
  &conf_fac(i,j+2,k+3) +    conf_fac(i,j+3,k+3) &
                    +  9*conf_fac(i,j-3,k+2) -  81*conf_fac(i,j-2,k+2) +  405*conf_fac(i,j-1,k+2) -  405*conf_fac(i,j+1,k+2) +  81*&
  &conf_fac(i,j+2,k+2) -  9*conf_fac(i,j+3,k+2) &
                    - 45*conf_fac(i,j-3,k+1) + 405*conf_fac(i,j-2,k+1) - 2025*conf_fac(i,j-1,k+1) + 2025*conf_fac(i,j+1,k+1) - 405*&
  &conf_fac(i,j+2,k+1) + 45*conf_fac(i,j+3,k+1) &
                    + 45*conf_fac(i,j-3,k-1) - 405*conf_fac(i,j-2,k-1) + 2025*conf_fac(i,j-1,k-1) - 2025*conf_fac(i,j+1,k-1) + 405*&
  &conf_fac(i,j+2,k-1) - 45*conf_fac(i,j+3,k-1) &
                    -  9*conf_fac(i,j-3,k-2) +  81*conf_fac(i,j-2,k-2) -  405*conf_fac(i,j-1,k-2) +  405*conf_fac(i,j+1,k-2) -  81*&
  &conf_fac(i,j+2,k-2) +  9*conf_fac(i,j+3,k-2) &
                    +    conf_fac(i,j-3,k-3) -   9*conf_fac(i,j-2,k-3) +   45*conf_fac(i,j-1,k-3) -   45*conf_fac(i,j+1,k-3) +   9*&
  &conf_fac(i,j+2,k-3) -    conf_fac(i,j+3,k-3) ) * odydz3600

      d2_ww(2,1) = d2_ww(1,2)
      d2_ww(3,1) = d2_ww(1,3)
      d2_ww(3,2) = d2_ww(2,3)


      
      d2_hh11(1,1) = (  2*hxx(i+3,j,k) - 27*hxx(i+2,j,k) + 270*hxx(i+1,j,k) - 490*hxx(i,j,k)&
                      + 2*hxx(i-3,j,k) - 27*hxx(i-2,j,k) + 270*hxx(i-1,j,k) ) * odxsq180
      d2_hh12(1,1) = (  2*hxy(i+3,j,k) - 27*hxy(i+2,j,k) + 270*hxy(i+1,j,k) - 490*hxy(i,j,k)&
                      + 2*hxy(i-3,j,k) - 27*hxy(i-2,j,k) + 270*hxy(i-1,j,k) ) * odxsq180
      d2_hh13(1,1) = (  2*hxz(i+3,j,k) - 27*hxz(i+2,j,k) + 270*hxz(i+1,j,k) - 490*hxz(i,j,k)&
                      + 2*hxz(i-3,j,k) - 27*hxz(i-2,j,k) + 270*hxz(i-1,j,k) ) * odxsq180
      d2_hh22(1,1) = (  2*hyy(i+3,j,k) - 27*hyy(i+2,j,k) + 270*hyy(i+1,j,k) - 490*hyy(i,j,k)&
                      + 2*hyy(i-3,j,k) - 27*hyy(i-2,j,k) + 270*hyy(i-1,j,k) ) * odxsq180
      d2_hh23(1,1) = (  2*hyz(i+3,j,k) - 27*hyz(i+2,j,k) + 270*hyz(i+1,j,k) - 490*hyz(i,j,k)&
                      + 2*hyz(i-3,j,k) - 27*hyz(i-2,j,k) + 270*hyz(i-1,j,k) ) * odxsq180
      d2_hh33(1,1) = (  2*hzz(i+3,j,k) - 27*hzz(i+2,j,k) + 270*hzz(i+1,j,k) - 490*hzz(i,j,k)&
                      + 2*hzz(i-3,j,k) - 27*hzz(i-2,j,k) + 270*hzz(i-1,j,k) ) * odxsq180

      d2_hh11(2,2) = (  2*hxx(i,j+3,k) - 27*hxx(i,j+2,k) + 270*hxx(i,j+1,k) - 490*hxx(i,j,k)&
                      + 2*hxx(i,j-3,k) - 27*hxx(i,j-2,k) + 270*hxx(i,j-1,k) ) * odysq180
      d2_hh12(2,2) = (  2*hxy(i,j+3,k) - 27*hxy(i,j+2,k) + 270*hxy(i,j+1,k) - 490*hxy(i,j,k)&
                      + 2*hxy(i,j-3,k) - 27*hxy(i,j-2,k) + 270*hxy(i,j-1,k) ) * odysq180
      d2_hh13(2,2) = (  2*hxz(i,j+3,k) - 27*hxz(i,j+2,k) + 270*hxz(i,j+1,k) - 490*hxz(i,j,k)&
                      + 2*hxz(i,j-3,k) - 27*hxz(i,j-2,k) + 270*hxz(i,j-1,k) ) * odysq180
      d2_hh22(2,2) = (  2*hyy(i,j+3,k) - 27*hyy(i,j+2,k) + 270*hyy(i,j+1,k) - 490*hyy(i,j,k)&
                      + 2*hyy(i,j-3,k) - 27*hyy(i,j-2,k) + 270*hyy(i,j-1,k) ) * odysq180
      d2_hh23(2,2) = (  2*hyz(i,j+3,k) - 27*hyz(i,j+2,k) + 270*hyz(i,j+1,k) - 490*hyz(i,j,k)&
                      + 2*hyz(i,j-3,k) - 27*hyz(i,j-2,k) + 270*hyz(i,j-1,k) ) * odysq180
      d2_hh33(2,2) = (  2*hzz(i,j+3,k) - 27*hzz(i,j+2,k) + 270*hzz(i,j+1,k) - 490*hzz(i,j,k)&
                      + 2*hzz(i,j-3,k) - 27*hzz(i,j-2,k) + 270*hzz(i,j-1,k) ) * odysq180

      d2_hh11(3,3) = (  2*hxx(i,j,k+3) - 27*hxx(i,j,k+2) + 270*hxx(i,j,k+1) - 490*hxx(i,j,k)&
                      + 2*hxx(i,j,k-3) - 27*hxx(i,j,k-2) + 270*hxx(i,j,k-1) ) * odzsq180
      d2_hh12(3,3) = (  2*hxy(i,j,k+3) - 27*hxy(i,j,k+2) + 270*hxy(i,j,k+1) - 490*hxy(i,j,k)&
                      + 2*hxy(i,j,k-3) - 27*hxy(i,j,k-2) + 270*hxy(i,j,k-1) ) * odzsq180
      d2_hh13(3,3) = (  2*hxz(i,j,k+3) - 27*hxz(i,j,k+2) + 270*hxz(i,j,k+1) - 490*hxz(i,j,k)&
                      + 2*hxz(i,j,k-3) - 27*hxz(i,j,k-2) + 270*hxz(i,j,k-1) ) * odzsq180
      d2_hh22(3,3) = (  2*hyy(i,j,k+3) - 27*hyy(i,j,k+2) + 270*hyy(i,j,k+1) - 490*hyy(i,j,k)&
                      + 2*hyy(i,j,k-3) - 27*hyy(i,j,k-2) + 270*hyy(i,j,k-1) ) * odzsq180
      d2_hh23(3,3) = (  2*hyz(i,j,k+3) - 27*hyz(i,j,k+2) + 270*hyz(i,j,k+1) - 490*hyz(i,j,k)&
                      + 2*hyz(i,j,k-3) - 27*hyz(i,j,k-2) + 270*hyz(i,j,k-1) ) * odzsq180
      d2_hh33(3,3) = (  2*hzz(i,j,k+3) - 27*hzz(i,j,k+2) + 270*hzz(i,j,k+1) - 490*hzz(i,j,k)&
                      + 2*hzz(i,j,k-3) - 27*hzz(i,j,k-2) + 270*hzz(i,j,k-1) ) * odzsq180

      d2_hh11(1,2) = (    -hxx(i-3,j+3,k) +   9*hxx(i-2,j+3,k) -   45*hxx(i-1,j+3,k) +   45*hxx(i+1,j+3,k) -   9*hxx(i+2,j+3,k) +  &
  &  hxx(i+3,j+3,k) &
                      +  9*hxx(i-3,j+2,k) -  81*hxx(i-2,j+2,k) +  405*hxx(i-1,j+2,k) -  405*hxx(i+1,j+2,k) +  81*hxx(i+2,j+2,k) -  &
  &9*hxx(i+3,j+2,k) &
                      - 45*hxx(i-3,j+1,k) + 405*hxx(i-2,j+1,k) - 2025*hxx(i-1,j+1,k) + 2025*hxx(i+1,j+1,k) - 405*hxx(i+2,j+1,k) + 4&
  &5*hxx(i+3,j+1,k) &
                      + 45*hxx(i-3,j-1,k) - 405*hxx(i-2,j-1,k) + 2025*hxx(i-1,j-1,k) - 2025*hxx(i+1,j-1,k) + 405*hxx(i+2,j-1,k) - 4&
  &5*hxx(i+3,j-1,k) &
                      -  9*hxx(i-3,j-2,k) +  81*hxx(i-2,j-2,k) -  405*hxx(i-1,j-2,k) +  405*hxx(i+1,j-2,k) -  81*hxx(i+2,j-2,k) +  &
  &9*hxx(i+3,j-2,k) &
                      +    hxx(i-3,j-3,k) -   9*hxx(i-2,j-3,k) +   45*hxx(i-1,j-3,k) -   45*hxx(i+1,j-3,k) +   9*hxx(i+2,j-3,k) -  &
  &  hxx(i+3,j-3,k) ) * odxdy3600
      d2_hh12(1,2) = (    -hxy(i-3,j+3,k) +   9*hxy(i-2,j+3,k) -   45*hxy(i-1,j+3,k) +   45*hxy(i+1,j+3,k) -   9*hxy(i+2,j+3,k) +  &
  &  hxy(i+3,j+3,k) &
                      +  9*hxy(i-3,j+2,k) -  81*hxy(i-2,j+2,k) +  405*hxy(i-1,j+2,k) -  405*hxy(i+1,j+2,k) +  81*hxy(i+2,j+2,k) -  &
  &9*hxy(i+3,j+2,k) &
                      - 45*hxy(i-3,j+1,k) + 405*hxy(i-2,j+1,k) - 2025*hxy(i-1,j+1,k) + 2025*hxy(i+1,j+1,k) - 405*hxy(i+2,j+1,k) + 4&
  &5*hxy(i+3,j+1,k) &
                      + 45*hxy(i-3,j-1,k) - 405*hxy(i-2,j-1,k) + 2025*hxy(i-1,j-1,k) - 2025*hxy(i+1,j-1,k) + 405*hxy(i+2,j-1,k) - 4&
  &5*hxy(i+3,j-1,k) &
                      -  9*hxy(i-3,j-2,k) +  81*hxy(i-2,j-2,k) -  405*hxy(i-1,j-2,k) +  405*hxy(i+1,j-2,k) -  81*hxy(i+2,j-2,k) +  &
  &9*hxy(i+3,j-2,k) &
                      +    hxy(i-3,j-3,k) -   9*hxy(i-2,j-3,k) +   45*hxy(i-1,j-3,k) -   45*hxy(i+1,j-3,k) +   9*hxy(i+2,j-3,k) -  &
  &  hxy(i+3,j-3,k) ) * odxdy3600
      d2_hh13(1,2) = (    -hxz(i-3,j+3,k) +   9*hxz(i-2,j+3,k) -   45*hxz(i-1,j+3,k) +   45*hxz(i+1,j+3,k) -   9*hxz(i+2,j+3,k) +  &
  &  hxz(i+3,j+3,k) &
                      +  9*hxz(i-3,j+2,k) -  81*hxz(i-2,j+2,k) +  405*hxz(i-1,j+2,k) -  405*hxz(i+1,j+2,k) +  81*hxz(i+2,j+2,k) -  &
  &9*hxz(i+3,j+2,k) &
                      - 45*hxz(i-3,j+1,k) + 405*hxz(i-2,j+1,k) - 2025*hxz(i-1,j+1,k) + 2025*hxz(i+1,j+1,k) - 405*hxz(i+2,j+1,k) + 4&
  &5*hxz(i+3,j+1,k) &
                      + 45*hxz(i-3,j-1,k) - 405*hxz(i-2,j-1,k) + 2025*hxz(i-1,j-1,k) - 2025*hxz(i+1,j-1,k) + 405*hxz(i+2,j-1,k) - 4&
  &5*hxz(i+3,j-1,k) &
                      -  9*hxz(i-3,j-2,k) +  81*hxz(i-2,j-2,k) -  405*hxz(i-1,j-2,k) +  405*hxz(i+1,j-2,k) -  81*hxz(i+2,j-2,k) +  &
  &9*hxz(i+3,j-2,k) &
                      +    hxz(i-3,j-3,k) -   9*hxz(i-2,j-3,k) +   45*hxz(i-1,j-3,k) -   45*hxz(i+1,j-3,k) +   9*hxz(i+2,j-3,k) -  &
  &  hxz(i+3,j-3,k) ) * odxdy3600
      d2_hh22(1,2) = (    -hyy(i-3,j+3,k) +   9*hyy(i-2,j+3,k) -   45*hyy(i-1,j+3,k) +   45*hyy(i+1,j+3,k) -   9*hyy(i+2,j+3,k) +  &
  &  hyy(i+3,j+3,k) &
                      +  9*hyy(i-3,j+2,k) -  81*hyy(i-2,j+2,k) +  405*hyy(i-1,j+2,k) -  405*hyy(i+1,j+2,k) +  81*hyy(i+2,j+2,k) -  &
  &9*hyy(i+3,j+2,k) &
                      - 45*hyy(i-3,j+1,k) + 405*hyy(i-2,j+1,k) - 2025*hyy(i-1,j+1,k) + 2025*hyy(i+1,j+1,k) - 405*hyy(i+2,j+1,k) + 4&
  &5*hyy(i+3,j+1,k) &
                      + 45*hyy(i-3,j-1,k) - 405*hyy(i-2,j-1,k) + 2025*hyy(i-1,j-1,k) - 2025*hyy(i+1,j-1,k) + 405*hyy(i+2,j-1,k) - 4&
  &5*hyy(i+3,j-1,k) &
                      -  9*hyy(i-3,j-2,k) +  81*hyy(i-2,j-2,k) -  405*hyy(i-1,j-2,k) +  405*hyy(i+1,j-2,k) -  81*hyy(i+2,j-2,k) +  &
  &9*hyy(i+3,j-2,k) &
                      +    hyy(i-3,j-3,k) -   9*hyy(i-2,j-3,k) +   45*hyy(i-1,j-3,k) -   45*hyy(i+1,j-3,k) +   9*hyy(i+2,j-3,k) -  &
  &  hyy(i+3,j-3,k) ) * odxdy3600
      d2_hh23(1,2) = (    -hyz(i-3,j+3,k) +   9*hyz(i-2,j+3,k) -   45*hyz(i-1,j+3,k) +   45*hyz(i+1,j+3,k) -   9*hyz(i+2,j+3,k) +  &
  &  hyz(i+3,j+3,k) &
                      +  9*hyz(i-3,j+2,k) -  81*hyz(i-2,j+2,k) +  405*hyz(i-1,j+2,k) -  405*hyz(i+1,j+2,k) +  81*hyz(i+2,j+2,k) -  &
  &9*hyz(i+3,j+2,k) &
                      - 45*hyz(i-3,j+1,k) + 405*hyz(i-2,j+1,k) - 2025*hyz(i-1,j+1,k) + 2025*hyz(i+1,j+1,k) - 405*hyz(i+2,j+1,k) + 4&
  &5*hyz(i+3,j+1,k) &
                      + 45*hyz(i-3,j-1,k) - 405*hyz(i-2,j-1,k) + 2025*hyz(i-1,j-1,k) - 2025*hyz(i+1,j-1,k) + 405*hyz(i+2,j-1,k) - 4&
  &5*hyz(i+3,j-1,k) &
                      -  9*hyz(i-3,j-2,k) +  81*hyz(i-2,j-2,k) -  405*hyz(i-1,j-2,k) +  405*hyz(i+1,j-2,k) -  81*hyz(i+2,j-2,k) +  &
  &9*hyz(i+3,j-2,k) &
                      +    hyz(i-3,j-3,k) -   9*hyz(i-2,j-3,k) +   45*hyz(i-1,j-3,k) -   45*hyz(i+1,j-3,k) +   9*hyz(i+2,j-3,k) -  &
  &  hyz(i+3,j-3,k) ) * odxdy3600
      d2_hh33(1,2) = (    -hzz(i-3,j+3,k) +   9*hzz(i-2,j+3,k) -   45*hzz(i-1,j+3,k) +   45*hzz(i+1,j+3,k) -   9*hzz(i+2,j+3,k) +  &
  &  hzz(i+3,j+3,k) &
                      +  9*hzz(i-3,j+2,k) -  81*hzz(i-2,j+2,k) +  405*hzz(i-1,j+2,k) -  405*hzz(i+1,j+2,k) +  81*hzz(i+2,j+2,k) -  &
  &9*hzz(i+3,j+2,k) &
                      - 45*hzz(i-3,j+1,k) + 405*hzz(i-2,j+1,k) - 2025*hzz(i-1,j+1,k) + 2025*hzz(i+1,j+1,k) - 405*hzz(i+2,j+1,k) + 4&
  &5*hzz(i+3,j+1,k) &
                      + 45*hzz(i-3,j-1,k) - 405*hzz(i-2,j-1,k) + 2025*hzz(i-1,j-1,k) - 2025*hzz(i+1,j-1,k) + 405*hzz(i+2,j-1,k) - 4&
  &5*hzz(i+3,j-1,k) &
                      -  9*hzz(i-3,j-2,k) +  81*hzz(i-2,j-2,k) -  405*hzz(i-1,j-2,k) +  405*hzz(i+1,j-2,k) -  81*hzz(i+2,j-2,k) +  &
  &9*hzz(i+3,j-2,k) &
                      +    hzz(i-3,j-3,k) -   9*hzz(i-2,j-3,k) +   45*hzz(i-1,j-3,k) -   45*hzz(i+1,j-3,k) +   9*hzz(i+2,j-3,k) -  &
  &  hzz(i+3,j-3,k) ) * odxdy3600

      d2_hh11(1,3) = (    -hxx(i-3,j,k+3) +   9*hxx(i-2,j,k+3) -   45*hxx(i-1,j,k+3) +   45*hxx(i+1,j,k+3) -   9*hxx(i+2,j,k+3) +  &
  &  hxx(i+3,j,k+3) &
                      +  9*hxx(i-3,j,k+2) -  81*hxx(i-2,j,k+2) +  405*hxx(i-1,j,k+2) -  405*hxx(i+1,j,k+2) +  81*hxx(i+2,j,k+2) -  &
  &9*hxx(i+3,j,k+2) &
                      - 45*hxx(i-3,j,k+1) + 405*hxx(i-2,j,k+1) - 2025*hxx(i-1,j,k+1) + 2025*hxx(i+1,j,k+1) - 405*hxx(i+2,j,k+1) + 4&
  &5*hxx(i+3,j,k+1) &
                      + 45*hxx(i-3,j,k-1) - 405*hxx(i-2,j,k-1) + 2025*hxx(i-1,j,k-1) - 2025*hxx(i+1,j,k-1) + 405*hxx(i+2,j,k-1) - 4&
  &5*hxx(i+3,j,k-1) &
                      -  9*hxx(i-3,j,k-2) +  81*hxx(i-2,j,k-2) -  405*hxx(i-1,j,k-2) +  405*hxx(i+1,j,k-2) -  81*hxx(i+2,j,k-2) +  &
  &9*hxx(i+3,j,k-2) &
                      +    hxx(i-3,j,k-3) -   9*hxx(i-2,j,k-3) +   45*hxx(i-1,j,k-3) -   45*hxx(i+1,j,k-3) +   9*hxx(i+2,j,k-3) -  &
  &  hxx(i+3,j,k-3) ) * odxdz3600
      d2_hh12(1,3) = (    -hxy(i-3,j,k+3) +   9*hxy(i-2,j,k+3) -   45*hxy(i-1,j,k+3) +   45*hxy(i+1,j,k+3) -   9*hxy(i+2,j,k+3) +  &
  &  hxy(i+3,j,k+3) &
                      +  9*hxy(i-3,j,k+2) -  81*hxy(i-2,j,k+2) +  405*hxy(i-1,j,k+2) -  405*hxy(i+1,j,k+2) +  81*hxy(i+2,j,k+2) -  &
  &9*hxy(i+3,j,k+2) &
                        - 45*hxy(i-3,j,k+1) + 405*hxy(i-2,j,k+1) - 2025*hxy(i-1,j,k+1) + 2025*hxy(i+1,j,k+1) - 405*hxy(i+2,j,k+1) +&
  & 45*hxy(i+3,j,k+1) &
                      + 45*hxy(i-3,j,k-1) - 405*hxy(i-2,j,k-1) + 2025*hxy(i-1,j,k-1) - 2025*hxy(i+1,j,k-1) + 405*hxy(i+2,j,k-1) - 4&
  &5*hxy(i+3,j,k-1) &
                      -  9*hxy(i-3,j,k-2) +  81*hxy(i-2,j,k-2) -  405*hxy(i-1,j,k-2) +  405*hxy(i+1,j,k-2) -  81*hxy(i+2,j,k-2) +  &
  &9*hxy(i+3,j,k-2) &
                      +    hxy(i-3,j,k-3) -   9*hxy(i-2,j,k-3) +   45*hxy(i-1,j,k-3) -   45*hxy(i+1,j,k-3) +   9*hxy(i+2,j,k-3) -  &
  &  hxy(i+3,j,k-3) ) * odxdz3600
      d2_hh13(1,3) = (    -hxz(i-3,j,k+3) +   9*hxz(i-2,j,k+3) -   45*hxz(i-1,j,k+3) +   45*hxz(i+1,j,k+3) -   9*hxz(i+2,j,k+3) +  &
  &  hxz(i+3,j,k+3) &
                      +  9*hxz(i-3,j,k+2) -  81*hxz(i-2,j,k+2) +  405*hxz(i-1,j,k+2) -  405*hxz(i+1,j,k+2) +  81*hxz(i+2,j,k+2) -  &
  &9*hxz(i+3,j,k+2) &
                      - 45*hxz(i-3,j,k+1) + 405*hxz(i-2,j,k+1) - 2025*hxz(i-1,j,k+1) + 2025*hxz(i+1,j,k+1) - 405*hxz(i+2,j,k+1) + 4&
  &5*hxz(i+3,j,k+1) &
                      + 45*hxz(i-3,j,k-1) - 405*hxz(i-2,j,k-1) + 2025*hxz(i-1,j,k-1) - 2025*hxz(i+1,j,k-1) + 405*hxz(i+2,j,k-1) - 4&
  &5*hxz(i+3,j,k-1) &
                      -  9*hxz(i-3,j,k-2) +  81*hxz(i-2,j,k-2) -  405*hxz(i-1,j,k-2) +  405*hxz(i+1,j,k-2) -  81*hxz(i+2,j,k-2) +  &
  &9*hxz(i+3,j,k-2) &
                      +    hxz(i-3,j,k-3) -   9*hxz(i-2,j,k-3) +   45*hxz(i-1,j,k-3) -   45*hxz(i+1,j,k-3) +   9*hxz(i+2,j,k-3) -  &
  &  hxz(i+3,j,k-3) ) * odxdz3600
      d2_hh22(1,3) = (    -hyy(i-3,j,k+3) +   9*hyy(i-2,j,k+3) -   45*hyy(i-1,j,k+3) +   45*hyy(i+1,j,k+3) -   9*hyy(i+2,j,k+3) +  &
  &  hyy(i+3,j,k+3) &
                      +  9*hyy(i-3,j,k+2) -  81*hyy(i-2,j,k+2) +  405*hyy(i-1,j,k+2) -  405*hyy(i+1,j,k+2) +  81*hyy(i+2,j,k+2) -  &
  &9*hyy(i+3,j,k+2) &
                      - 45*hyy(i-3,j,k+1) + 405*hyy(i-2,j,k+1) - 2025*hyy(i-1,j,k+1) + 2025*hyy(i+1,j,k+1) - 405*hyy(i+2,j,k+1) + 4&
  &5*hyy(i+3,j,k+1) &
                      + 45*hyy(i-3,j,k-1) - 405*hyy(i-2,j,k-1) + 2025*hyy(i-1,j,k-1) - 2025*hyy(i+1,j,k-1) + 405*hyy(i+2,j,k-1) - 4&
  &5*hyy(i+3,j,k-1) &
                      -  9*hyy(i-3,j,k-2) +  81*hyy(i-2,j,k-2) -  405*hyy(i-1,j,k-2) +  405*hyy(i+1,j,k-2) -  81*hyy(i+2,j,k-2) +  &
  &9*hyy(i+3,j,k-2) &
                      +    hyy(i-3,j,k-3) -   9*hyy(i-2,j,k-3) +   45*hyy(i-1,j,k-3) -   45*hyy(i+1,j,k-3) +   9*hyy(i+2,j,k-3) -  &
  &  hyy(i+3,j,k-3) ) * odxdz3600
      d2_hh23(1,3) = (    -hyz(i-3,j,k+3) +   9*hyz(i-2,j,k+3) -   45*hyz(i-1,j,k+3) +   45*hyz(i+1,j,k+3) -   9*hyz(i+2,j,k+3) +  &
  &  hyz(i+3,j,k+3) &
                      +  9*hyz(i-3,j,k+2) -  81*hyz(i-2,j,k+2) +  405*hyz(i-1,j,k+2) -  405*hyz(i+1,j,k+2) +  81*hyz(i+2,j,k+2) -  &
  &9*hyz(i+3,j,k+2) &
                      - 45*hyz(i-3,j,k+1) + 405*hyz(i-2,j,k+1) - 2025*hyz(i-1,j,k+1) + 2025*hyz(i+1,j,k+1) - 405*hyz(i+2,j,k+1) + 4&
  &5*hyz(i+3,j,k+1) &
                      + 45*hyz(i-3,j,k-1) - 405*hyz(i-2,j,k-1) + 2025*hyz(i-1,j,k-1) - 2025*hyz(i+1,j,k-1) + 405*hyz(i+2,j,k-1) - 4&
  &5*hyz(i+3,j,k-1) &
                      -  9*hyz(i-3,j,k-2) +  81*hyz(i-2,j,k-2) -  405*hyz(i-1,j,k-2) +  405*hyz(i+1,j,k-2) -  81*hyz(i+2,j,k-2) +  &
  &9*hyz(i+3,j,k-2) &
                      +    hyz(i-3,j,k-3) -   9*hyz(i-2,j,k-3) +   45*hyz(i-1,j,k-3) -   45*hyz(i+1,j,k-3) +   9*hyz(i+2,j,k-3) -  &
  &  hyz(i+3,j,k-3) ) * odxdz3600
      d2_hh33(1,3) = (    -hzz(i-3,j,k+3) +   9*hzz(i-2,j,k+3) -   45*hzz(i-1,j,k+3) +   45*hzz(i+1,j,k+3) -   9*hzz(i+2,j,k+3) +  &
  &  hzz(i+3,j,k+3) &
                      +  9*hzz(i-3,j,k+2) -  81*hzz(i-2,j,k+2) +  405*hzz(i-1,j,k+2) -  405*hzz(i+1,j,k+2) +  81*hzz(i+2,j,k+2) -  &
  &9*hzz(i+3,j,k+2) &
                      - 45*hzz(i-3,j,k+1) + 405*hzz(i-2,j,k+1) - 2025*hzz(i-1,j,k+1) + 2025*hzz(i+1,j,k+1) - 405*hzz(i+2,j,k+1) + 4&
  &5*hzz(i+3,j,k+1) &
                      + 45*hzz(i-3,j,k-1) - 405*hzz(i-2,j,k-1) + 2025*hzz(i-1,j,k-1) - 2025*hzz(i+1,j,k-1) + 405*hzz(i+2,j,k-1) - 4&
  &5*hzz(i+3,j,k-1) &
                      -  9*hzz(i-3,j,k-2) +  81*hzz(i-2,j,k-2) -  405*hzz(i-1,j,k-2) +  405*hzz(i+1,j,k-2) -  81*hzz(i+2,j,k-2) +  &
  &9*hzz(i+3,j,k-2) &
                      +    hzz(i-3,j,k-3) -   9*hzz(i-2,j,k-3) +   45*hzz(i-1,j,k-3) -   45*hzz(i+1,j,k-3) +   9*hzz(i+2,j,k-3) -  &
  &  hzz(i+3,j,k-3) ) * odxdz3600

      d2_hh11(2,3) = (    -hxx(i,j-3,k+3) +   9*hxx(i,j-2,k+3) -   45*hxx(i,j-1,k+3) +   45*hxx(i,j+1,k+3) -   9*hxx(i,j+2,k+3) +  &
  &  hxx(i,j+3,k+3) &
                      +  9*hxx(i,j-3,k+2) -  81*hxx(i,j-2,k+2) +  405*hxx(i,j-1,k+2) -  405*hxx(i,j+1,k+2) +  81*hxx(i,j+2,k+2) -  &
  &9*hxx(i,j+3,k+2) &
                      - 45*hxx(i,j-3,k+1) + 405*hxx(i,j-2,k+1) - 2025*hxx(i,j-1,k+1) + 2025*hxx(i,j+1,k+1) - 405*hxx(i,j+2,k+1) + 4&
  &5*hxx(i,j+3,k+1) &
                      + 45*hxx(i,j-3,k-1) - 405*hxx(i,j-2,k-1) + 2025*hxx(i,j-1,k-1) - 2025*hxx(i,j+1,k-1) + 405*hxx(i,j+2,k-1) - 4&
  &5*hxx(i,j+3,k-1) &
                      -  9*hxx(i,j-3,k-2) +  81*hxx(i,j-2,k-2) -  405*hxx(i,j-1,k-2) +  405*hxx(i,j+1,k-2) -  81*hxx(i,j+2,k-2) +  &
  &9*hxx(i,j+3,k-2) &
                      +    hxx(i,j-3,k-3) -   9*hxx(i,j-2,k-3) +   45*hxx(i,j-1,k-3) -   45*hxx(i,j+1,k-3) +   9*hxx(i,j+2,k-3) -  &
  &  hxx(i,j+3,k-3) ) * odydz3600
      d2_hh12(2,3) = (    -hxy(i,j-3,k+3) +   9*hxy(i,j-2,k+3) -   45*hxy(i,j-1,k+3) +   45*hxy(i,j+1,k+3) -   9*hxy(i,j+2,k+3) +  &
  &  hxy(i,j+3,k+3) &
                      +  9*hxy(i,j-3,k+2) -  81*hxy(i,j-2,k+2) +  405*hxy(i,j-1,k+2) -  405*hxy(i,j+1,k+2) +  81*hxy(i,j+2,k+2) -  &
  &9*hxy(i,j+3,k+2) &
                      - 45*hxy(i,j-3,k+1) + 405*hxy(i,j-2,k+1) - 2025*hxy(i,j-1,k+1) + 2025*hxy(i,j+1,k+1) - 405*hxy(i,j+2,k+1) + 4&
  &5*hxy(i,j+3,k+1) &
                      + 45*hxy(i,j-3,k-1) - 405*hxy(i,j-2,k-1) + 2025*hxy(i,j-1,k-1) - 2025*hxy(i,j+1,k-1) + 405*hxy(i,j+2,k-1) - 4&
  &5*hxy(i,j+3,k-1) &
                      -  9*hxy(i,j-3,k-2) +  81*hxy(i,j-2,k-2) -  405*hxy(i,j-1,k-2) +  405*hxy(i,j+1,k-2) -  81*hxy(i,j+2,k-2) +  &
  &9*hxy(i,j+3,k-2) &
                      +    hxy(i,j-3,k-3) -   9*hxy(i,j-2,k-3) +   45*hxy(i,j-1,k-3) -   45*hxy(i,j+1,k-3) +   9*hxy(i,j+2,k-3) -  &
  &  hxy(i,j+3,k-3) ) * odydz3600
      d2_hh13(2,3) = (    -hxz(i,j-3,k+3) +   9*hxz(i,j-2,k+3) -   45*hxz(i,j-1,k+3) +   45*hxz(i,j+1,k+3) -   9*hxz(i,j+2,k+3) +  &
  &  hxz(i,j+3,k+3) &
                      +  9*hxz(i,j-3,k+2) -  81*hxz(i,j-2,k+2) +  405*hxz(i,j-1,k+2) -  405*hxz(i,j+1,k+2) +  81*hxz(i,j+2,k+2) -  &
  &9*hxz(i,j+3,k+2) &
                      - 45*hxz(i,j-3,k+1) + 405*hxz(i,j-2,k+1) - 2025*hxz(i,j-1,k+1) + 2025*hxz(i,j+1,k+1) - 405*hxz(i,j+2,k+1) + 4&
  &5*hxz(i,j+3,k+1) &
                      + 45*hxz(i,j-3,k-1) - 405*hxz(i,j-2,k-1) + 2025*hxz(i,j-1,k-1) - 2025*hxz(i,j+1,k-1) + 405*hxz(i,j+2,k-1) - 4&
  &5*hxz(i,j+3,k-1) &
                      -  9*hxz(i,j-3,k-2) +  81*hxz(i,j-2,k-2) -  405*hxz(i,j-1,k-2) +  405*hxz(i,j+1,k-2) -  81*hxz(i,j+2,k-2) +  &
  &9*hxz(i,j+3,k-2) &
                      +    hxz(i,j-3,k-3) -   9*hxz(i,j-2,k-3) +   45*hxz(i,j-1,k-3) -   45*hxz(i,j+1,k-3) +   9*hxz(i,j+2,k-3) -  &
  &  hxz(i,j+3,k-3) ) * odydz3600
      d2_hh22(2,3) = (    -hyy(i,j-3,k+3) +   9*hyy(i,j-2,k+3) -   45*hyy(i,j-1,k+3) +   45*hyy(i,j+1,k+3) -   9*hyy(i,j+2,k+3) +  &
  &  hyy(i,j+3,k+3) &
                      +  9*hyy(i,j-3,k+2) -  81*hyy(i,j-2,k+2) +  405*hyy(i,j-1,k+2) -  405*hyy(i,j+1,k+2) +  81*hyy(i,j+2,k+2) -  &
  &9*hyy(i,j+3,k+2) &
                      - 45*hyy(i,j-3,k+1) + 405*hyy(i,j-2,k+1) - 2025*hyy(i,j-1,k+1) + 2025*hyy(i,j+1,k+1) - 405*hyy(i,j+2,k+1) + 4&
  &5*hyy(i,j+3,k+1) &
                      + 45*hyy(i,j-3,k-1) - 405*hyy(i,j-2,k-1) + 2025*hyy(i,j-1,k-1) - 2025*hyy(i,j+1,k-1) + 405*hyy(i,j+2,k-1) - 4&
  &5*hyy(i,j+3,k-1) &
                      -  9*hyy(i,j-3,k-2) +  81*hyy(i,j-2,k-2) -  405*hyy(i,j-1,k-2) +  405*hyy(i,j+1,k-2) -  81*hyy(i,j+2,k-2) +  &
  &9*hyy(i,j+3,k-2) &
                      +    hyy(i,j-3,k-3) -   9*hyy(i,j-2,k-3) +   45*hyy(i,j-1,k-3) -   45*hyy(i,j+1,k-3) +   9*hyy(i,j+2,k-3) -  &
  &  hyy(i,j+3,k-3) ) * odydz3600
      d2_hh23(2,3) = (    -hyz(i,j-3,k+3) +   9*hyz(i,j-2,k+3) -   45*hyz(i,j-1,k+3) +   45*hyz(i,j+1,k+3) -   9*hyz(i,j+2,k+3) +  &
  &  hyz(i,j+3,k+3) &
                      +  9*hyz(i,j-3,k+2) -  81*hyz(i,j-2,k+2) +  405*hyz(i,j-1,k+2) -  405*hyz(i,j+1,k+2) +  81*hyz(i,j+2,k+2) -  &
  &9*hyz(i,j+3,k+2) &
                      - 45*hyz(i,j-3,k+1) + 405*hyz(i,j-2,k+1) - 2025*hyz(i,j-1,k+1) + 2025*hyz(i,j+1,k+1) - 405*hyz(i,j+2,k+1) + 4&
  &5*hyz(i,j+3,k+1) &
                      + 45*hyz(i,j-3,k-1) - 405*hyz(i,j-2,k-1) + 2025*hyz(i,j-1,k-1) - 2025*hyz(i,j+1,k-1) + 405*hyz(i,j+2,k-1) - 4&
  &5*hyz(i,j+3,k-1) &
                      -  9*hyz(i,j-3,k-2) +  81*hyz(i,j-2,k-2) -  405*hyz(i,j-1,k-2) +  405*hyz(i,j+1,k-2) -  81*hyz(i,j+2,k-2) +  &
  &9*hyz(i,j+3,k-2) &
                      +    hyz(i,j-3,k-3) -   9*hyz(i,j-2,k-3) +   45*hyz(i,j-1,k-3) -   45*hyz(i,j+1,k-3) +   9*hyz(i,j+2,k-3) -  &
  &  hyz(i,j+3,k-3) ) * odydz3600
      d2_hh33(2,3) = (    -hzz(i,j-3,k+3) +   9*hzz(i,j-2,k+3) -   45*hzz(i,j-1,k+3) +   45*hzz(i,j+1,k+3) -   9*hzz(i,j+2,k+3) +  &
  &  hzz(i,j+3,k+3) &
                        +  9*hzz(i,j-3,k+2) -  81*hzz(i,j-2,k+2) +  405*hzz(i,j-1,k+2) -  405*hzz(i,j+1,k+2) +  81*hzz(i,j+2,k+2) -&
  &  9*hzz(i,j+3,k+2) &
                        - 45*hzz(i,j-3,k+1) + 405*hzz(i,j-2,k+1) - 2025*hzz(i,j-1,k+1) + 2025*hzz(i,j+1,k+1) - 405*hzz(i,j+2,k+1) +&
  & 45*hzz(i,j+3,k+1) &
                        + 45*hzz(i,j-3,k-1) - 405*hzz(i,j-2,k-1) + 2025*hzz(i,j-1,k-1) - 2025*hzz(i,j+1,k-1) + 405*hzz(i,j+2,k-1) -&
  & 45*hzz(i,j+3,k-1) &
                        -  9*hzz(i,j-3,k-2) +  81*hzz(i,j-2,k-2) -  405*hzz(i,j-1,k-2) +  405*hzz(i,j+1,k-2) -  81*hzz(i,j+2,k-2) +&
  &  9*hzz(i,j+3,k-2) &
                        +    hzz(i,j-3,k-3) -   9*hzz(i,j-2,k-3) +   45*hzz(i,j-1,k-3) -   45*hzz(i,j+1,k-3) +   9*hzz(i,j+2,k-3) -&
  &    hzz(i,j+3,k-3) ) * odydz3600

      d2_hh11(2,1) = d2_hh11(1,2)
      d2_hh12(2,1) = d2_hh12(1,2)
      d2_hh13(2,1) = d2_hh13(1,2)
      d2_hh22(2,1) = d2_hh22(1,2)
      d2_hh23(2,1) = d2_hh23(1,2)
      d2_hh33(2,1) = d2_hh33(1,2)

      d2_hh11(3,1) = d2_hh11(1,3)
      d2_hh12(3,1) = d2_hh12(1,3)
      d2_hh13(3,1) = d2_hh13(1,3)
      d2_hh22(3,1) = d2_hh22(1,3)
      d2_hh23(3,1) = d2_hh23(1,3)
      d2_hh33(3,1) = d2_hh33(1,3)

      d2_hh11(3,2) = d2_hh11(2,3)
      d2_hh12(3,2) = d2_hh12(2,3)
      d2_hh13(3,2) = d2_hh13(2,3)
      d2_hh22(3,2) = d2_hh22(2,3)
      d2_hh23(3,2) = d2_hh23(2,3)
      d2_hh33(3,2) = d2_hh33(2,3)

    else 

      
      
      d1_ww(1) = (   -conf_fac(i+2,j,k) + 8*conf_fac(i+1,j,k)                          &
                  - 8*conf_fac(i-1,j,k) +   conf_fac(i-2,j,k) ) / dx12

      d1_ww(2) = (   -conf_fac(i,j+2,k) + 8*conf_fac(i,j+1,k)                          &
                  - 8*conf_fac(i,j-1,k) +   conf_fac(i,j-2,k) ) / dy12

      d1_ww(3) = (   -conf_fac(i,j,k+2) + 8*conf_fac(i,j,k+1)                          &
                  - 8*conf_fac(i,j,k-1) +   conf_fac(i,j,k-2) ) / dz12

      
      d1_hh11(1) = (   -hxx(i+2,j,k) + 8*hxx(i+1,j,k)                      &
                    - 8*hxx(i-1,j,k) +   hxx(i-2,j,k) ) / dx12
      d1_hh12(1) = (   -hxy(i+2,j,k) + 8*hxy(i+1,j,k)                      &
                    - 8*hxy(i-1,j,k) +   hxy(i-2,j,k) ) / dx12
      d1_hh13(1) = (   -hxz(i+2,j,k) + 8*hxz(i+1,j,k)                      &
                    - 8*hxz(i-1,j,k) +   hxz(i-2,j,k) ) / dx12
      d1_hh22(1) = (   -hyy(i+2,j,k) + 8*hyy(i+1,j,k)                      &
                    - 8*hyy(i-1,j,k) +   hyy(i-2,j,k) ) / dx12
      d1_hh23(1) = (   -hyz(i+2,j,k) + 8*hyz(i+1,j,k)                      &
                    - 8*hyz(i-1,j,k) +   hyz(i-2,j,k) ) / dx12
      d1_hh33(1) = (   -hzz(i+2,j,k) + 8*hzz(i+1,j,k)                      &
                    - 8*hzz(i-1,j,k) +   hzz(i-2,j,k) ) / dx12

      d1_hh11(2) = (   -hxx(i,j+2,k) + 8*hxx(i,j+1,k)                      &
                    - 8*hxx(i,j-1,k) +   hxx(i,j-2,k) ) / dy12
      d1_hh12(2) = (   -hxy(i,j+2,k) + 8*hxy(i,j+1,k)                      &
                    - 8*hxy(i,j-1,k) +   hxy(i,j-2,k) ) / dy12
      d1_hh13(2) = (   -hxz(i,j+2,k) + 8*hxz(i,j+1,k)                      &
                    - 8*hxz(i,j-1,k) +   hxz(i,j-2,k) ) / dy12
      d1_hh22(2) = (   -hyy(i,j+2,k) + 8*hyy(i,j+1,k)                      &
                    - 8*hyy(i,j-1,k) +   hyy(i,j-2,k) ) / dy12
      d1_hh23(2) = (   -hyz(i,j+2,k) + 8*hyz(i,j+1,k)                      &
                    - 8*hyz(i,j-1,k) +   hyz(i,j-2,k) ) / dy12
      d1_hh33(2) = (   -hzz(i,j+2,k) + 8*hzz(i,j+1,k)                      &
                    - 8*hzz(i,j-1,k) +   hzz(i,j-2,k) ) / dy12

      d1_hh11(3) = (   -hxx(i,j,k+2) + 8*hxx(i,j,k+1)                      &
                    - 8*hxx(i,j,k-1) +   hxx(i,j,k-2) ) / dz12
      d1_hh12(3) = (   -hxy(i,j,k+2) + 8*hxy(i,j,k+1)                      &
                    - 8*hxy(i,j,k-1) +   hxy(i,j,k-2) ) / dz12
      d1_hh13(3) = (   -hxz(i,j,k+2) + 8*hxz(i,j,k+1)                      &
                    - 8*hxz(i,j,k-1) +   hxz(i,j,k-2) ) / dz12
      d1_hh22(3) = (   -hyy(i,j,k+2) + 8*hyy(i,j,k+1)                      &
                    - 8*hyy(i,j,k-1) +   hyy(i,j,k-2) ) / dz12
      d1_hh23(3) = (   -hyz(i,j,k+2) + 8*hyz(i,j,k+1)                      &
                    - 8*hyz(i,j,k-1) +   hyz(i,j,k-2) ) / dz12
      d1_hh33(3) = (   -hzz(i,j,k+2) + 8*hzz(i,j,k+1)                      &
                    - 8*hzz(i,j,k-1) +   hzz(i,j,k-2) ) / dz12

      
      d1_trk(1) = (   -tracek(i+2,j,k) + 8*tracek(i+1,j,k)                   &
                   - 8*tracek(i-1,j,k) +   tracek(i-2,j,k) ) / dx12

      d1_trk(2) = (   -tracek(i,j+2,k) + 8*tracek(i,j+1,k)                   &
                   - 8*tracek(i,j-1,k) +   tracek(i,j-2,k) ) / dy12

      d1_trk(3) = (   -tracek(i,j,k+2) + 8*tracek(i,j,k+1)                   &
                   - 8*tracek(i,j,k-1) +   tracek(i,j,k-2) ) / dz12

      
      d1_aa11(1) = (   -axx(i+2,j,k) + 8*axx(i+1,j,k)                      &
                    - 8*axx(i-1,j,k) +   axx(i-2,j,k) ) / dx12
      d1_aa12(1) = (   -axy(i+2,j,k) + 8*axy(i+1,j,k)                      &
                    - 8*axy(i-1,j,k) +   axy(i-2,j,k) ) / dx12
      d1_aa13(1) = (   -axz(i+2,j,k) + 8*axz(i+1,j,k)                      &
                    - 8*axz(i-1,j,k) +   axz(i-2,j,k) ) / dx12
      d1_aa22(1) = (   -ayy(i+2,j,k) + 8*ayy(i+1,j,k)                      &
                    - 8*ayy(i-1,j,k) +   ayy(i-2,j,k) ) / dx12
      d1_aa23(1) = (   -ayz(i+2,j,k) + 8*ayz(i+1,j,k)                      &
                    - 8*ayz(i-1,j,k) +   ayz(i-2,j,k) ) / dx12
      d1_aa33(1) = (   -azz(i+2,j,k) + 8*azz(i+1,j,k)                      &
                    - 8*azz(i-1,j,k) +   azz(i-2,j,k) ) / dx12

      d1_aa11(2) = (   -axx(i,j+2,k) + 8*axx(i,j+1,k)                      &
                    - 8*axx(i,j-1,k) +   axx(i,j-2,k) ) / dy12
      d1_aa12(2) = (   -axy(i,j+2,k) + 8*axy(i,j+1,k)                      &
                    - 8*axy(i,j-1,k) +   axy(i,j-2,k) ) / dy12
      d1_aa13(2) = (   -axz(i,j+2,k) + 8*axz(i,j+1,k)                      &
                    - 8*axz(i,j-1,k) +   axz(i,j-2,k) ) / dy12
      d1_aa22(2) = (   -ayy(i,j+2,k) + 8*ayy(i,j+1,k)                      &
                    - 8*ayy(i,j-1,k) +   ayy(i,j-2,k) ) / dy12
      d1_aa23(2) = (   -ayz(i,j+2,k) + 8*ayz(i,j+1,k)                      &
                    - 8*ayz(i,j-1,k) +   ayz(i,j-2,k) ) / dy12
      d1_aa33(2) = (   -azz(i,j+2,k) + 8*azz(i,j+1,k)                      &
                    - 8*azz(i,j-1,k) +   azz(i,j-2,k) ) / dy12

      d1_aa11(3) = (   -axx(i,j,k+2) + 8*axx(i,j,k+1)                      &
                    - 8*axx(i,j,k-1) +   axx(i,j,k-2) ) / dz12
      d1_aa12(3) = (   -axy(i,j,k+2) + 8*axy(i,j,k+1)                      &
                    - 8*axy(i,j,k-1) +   axy(i,j,k-2) ) / dz12
      d1_aa13(3) = (   -axz(i,j,k+2) + 8*axz(i,j,k+1)                      &
                    - 8*axz(i,j,k-1) +   axz(i,j,k-2) ) / dz12
      d1_aa22(3) = (   -ayy(i,j,k+2) + 8*ayy(i,j,k+1)                      &
                    - 8*ayy(i,j,k-1) +   ayy(i,j,k-2) ) / dz12
      d1_aa23(3) = (   -ayz(i,j,k+2) + 8*ayz(i,j,k+1)                      &
                    - 8*ayz(i,j,k-1) +   ayz(i,j,k-2) ) / dz12
      d1_aa33(3) = (   -azz(i,j,k+2) + 8*azz(i,j,k+1)                      &
                    - 8*azz(i,j,k-1) +   azz(i,j,k-2) ) / dz12

      
      d1_gammat1(1) = (   -gammatx(i+2,j,k) + 8*gammatx(i+1,j,k)               &
                       - 8*gammatx(i-1,j,k) +   gammatx(i-2,j,k) ) / dx12
      d1_gammat2(1) = (   -gammaty(i+2,j,k) + 8*gammaty(i+1,j,k)               &
                       - 8*gammaty(i-1,j,k) +   gammaty(i-2,j,k) ) / dx12
      d1_gammat3(1) = (   -gammatz(i+2,j,k) + 8*gammatz(i+1,j,k)               &
                       - 8*gammatz(i-1,j,k) +   gammatz(i-2,j,k) ) / dx12

      d1_gammat1(2) = (   -gammatx(i,j+2,k) + 8*gammatx(i,j+1,k)               &
                       - 8*gammatx(i,j-1,k) +   gammatx(i,j-2,k) ) / dy12
      d1_gammat2(2) = (   -gammaty(i,j+2,k) + 8*gammaty(i,j+1,k)               &
                       - 8*gammaty(i,j-1,k) +   gammaty(i,j-2,k) ) / dy12
      d1_gammat3(2) = (   -gammatz(i,j+2,k) + 8*gammatz(i,j+1,k)               &
                       - 8*gammatz(i,j-1,k) +   gammatz(i,j-2,k) ) / dy12

      d1_gammat1(3) = (   -gammatx(i,j,k+2) + 8*gammatx(i,j,k+1)               &
                       - 8*gammatx(i,j,k-1) +   gammatx(i,j,k-2) ) / dz12
      d1_gammat2(3) = (   -gammaty(i,j,k+2) + 8*gammaty(i,j,k+1)               &
                       - 8*gammaty(i,j,k-1) +   gammaty(i,j,k-2) ) / dz12
      d1_gammat3(3) = (   -gammatz(i,j,k+2) + 8*gammatz(i,j,k+1)               &
                       - 8*gammatz(i,j,k-1) +   gammatz(i,j,k-2) ) / dz12

      


      

      
      d2_ww(1,1) = (   -conf_fac(i+2,j,k) + 16*conf_fac(i+1,j,k) - 30*conf_fac(i,j,k)       &
                   + 16*conf_fac(i-1,j,k) -    conf_fac(i-2,j,k) ) / dxsq12

      d2_ww(2,2) = (   -conf_fac(i,j+2,k) + 16*conf_fac(i,j+1,k) - 30*conf_fac(i,j,k)       &
                   + 16*conf_fac(i,j-1,k) -    conf_fac(i,j-2,k) ) / dysq12

      d2_ww(3,3) = (   -conf_fac(i,j,k+2) + 16*conf_fac(i,j,k+1) - 30*conf_fac(i,j,k)       &
                   + 16*conf_fac(i,j,k-1) -    conf_fac(i,j,k-2) ) / dzsq12

      d2_ww(1,2) = (   -conf_fac(i-2,j+2,k) +  8*conf_fac(i-1,j+2,k) -  8*conf_fac(i+1,j+2,k) +   conf_fac(i+2,j+2,k)   &
                    + 8*conf_fac(i-2,j+1,k) - 64*conf_fac(i-1,j+1,k) + 64*conf_fac(i+1,j+1,k) - 8*conf_fac(i+2,j+1,k)   &
                    - 8*conf_fac(i-2,j-1,k) + 64*conf_fac(i-1,j-1,k) - 64*conf_fac(i+1,j-1,k) + 8*conf_fac(i+2,j-1,k)   &
                    +   conf_fac(i-2,j-2,k) -  8*conf_fac(i-1,j-2,k) +  8*conf_fac(i+1,j-2,k) -   conf_fac(i+2,j-2,k) ) / dxdy144

      d2_ww(1,3) = (   -conf_fac(i-2,j,k+2) +  8*conf_fac(i-1,j,k+2) -  8*conf_fac(i+1,j,k+2) +   conf_fac(i+2,j,k+2)   &
                    + 8*conf_fac(i-2,j,k+1) - 64*conf_fac(i-1,j,k+1) + 64*conf_fac(i+1,j,k+1) - 8*conf_fac(i+2,j,k+1)   &
                    - 8*conf_fac(i-2,j,k-1) + 64*conf_fac(i-1,j,k-1) - 64*conf_fac(i+1,j,k-1) + 8*conf_fac(i+2,j,k-1)   &
                    +   conf_fac(i-2,j,k-2) -  8*conf_fac(i-1,j,k-2) +  8*conf_fac(i+1,j,k-2) -   conf_fac(i+2,j,k-2) ) / dxdz144

      d2_ww(2,3) = (   -conf_fac(i,j-2,k+2) +  8*conf_fac(i,j-1,k+2) -  8*conf_fac(i,j+1,k+2) +   conf_fac(i,j+2,k+2)   &
                    + 8*conf_fac(i,j-2,k+1) - 64*conf_fac(i,j-1,k+1) + 64*conf_fac(i,j+1,k+1) - 8*conf_fac(i,j+2,k+1)   &
                    - 8*conf_fac(i,j-2,k-1) + 64*conf_fac(i,j-1,k-1) - 64*conf_fac(i,j+1,k-1) + 8*conf_fac(i,j+2,k-1)   &
                    +   conf_fac(i,j-2,k-2) -  8*conf_fac(i,j-1,k-2) +  8*conf_fac(i,j+1,k-2) -   conf_fac(i,j+2,k-2) ) / dydz144

      d2_ww(2,1) = d2_ww(1,2)
      d2_ww(3,1) = d2_ww(1,3)
      d2_ww(3,2) = d2_ww(2,3)

      
      d2_hh11(1,1) = (   -hxx(i+2,j,k) + 16*hxx(i+1,j,k) - 30*hxx(i,j,k)   &
                     + 16*hxx(i-1,j,k) -    hxx(i-2,j,k) ) / dxsq12
      d2_hh12(1,1) = (   -hxy(i+2,j,k) + 16*hxy(i+1,j,k) - 30*hxy(i,j,k)   &
                     + 16*hxy(i-1,j,k) -    hxy(i-2,j,k) ) / dxsq12
      d2_hh13(1,1) = (   -hxz(i+2,j,k) + 16*hxz(i+1,j,k) - 30*hxz(i,j,k)   &
                     + 16*hxz(i-1,j,k) -    hxz(i-2,j,k) ) / dxsq12
      d2_hh22(1,1) = (   -hyy(i+2,j,k) + 16*hyy(i+1,j,k) - 30*hyy(i,j,k)   &
                     + 16*hyy(i-1,j,k) -    hyy(i-2,j,k) ) / dxsq12
      d2_hh23(1,1) = (   -hyz(i+2,j,k) + 16*hyz(i+1,j,k) - 30*hyz(i,j,k)   &
                     + 16*hyz(i-1,j,k) -    hyz(i-2,j,k) ) / dxsq12
      d2_hh33(1,1) = (   -hzz(i+2,j,k) + 16*hzz(i+1,j,k) - 30*hzz(i,j,k)   &
                     + 16*hzz(i-1,j,k) -    hzz(i-2,j,k) ) / dxsq12

      d2_hh11(2,2) = (   -hxx(i,j+2,k) + 16*hxx(i,j+1,k) - 30*hxx(i,j,k)   &
                     + 16*hxx(i,j-1,k) -    hxx(i,j-2,k) ) / dysq12
      d2_hh12(2,2) = (   -hxy(i,j+2,k) + 16*hxy(i,j+1,k) - 30*hxy(i,j,k)   &
                     + 16*hxy(i,j-1,k) -    hxy(i,j-2,k) ) / dysq12
      d2_hh13(2,2) = (   -hxz(i,j+2,k) + 16*hxz(i,j+1,k) - 30*hxz(i,j,k)   &
                     + 16*hxz(i,j-1,k) -    hxz(i,j-2,k) ) / dysq12
      d2_hh22(2,2) = (   -hyy(i,j+2,k) + 16*hyy(i,j+1,k) - 30*hyy(i,j,k)   &
                     + 16*hyy(i,j-1,k) -    hyy(i,j-2,k) ) / dysq12
      d2_hh23(2,2) = (   -hyz(i,j+2,k) + 16*hyz(i,j+1,k) - 30*hyz(i,j,k)   &
                     + 16*hyz(i,j-1,k) -    hyz(i,j-2,k) ) / dysq12
      d2_hh33(2,2) = (   -hzz(i,j+2,k) + 16*hzz(i,j+1,k) - 30*hzz(i,j,k)   &
                     + 16*hzz(i,j-1,k) -    hzz(i,j-2,k) ) / dysq12

      d2_hh11(3,3) = (   -hxx(i,j,k+2) + 16*hxx(i,j,k+1) - 30*hxx(i,j,k)   &
                     + 16*hxx(i,j,k-1) -    hxx(i,j,k-2) ) / dzsq12
      d2_hh12(3,3) = (   -hxy(i,j,k+2) + 16*hxy(i,j,k+1) - 30*hxy(i,j,k)   &
                     + 16*hxy(i,j,k-1) -    hxy(i,j,k-2) ) / dzsq12
      d2_hh13(3,3) = (   -hxz(i,j,k+2) + 16*hxz(i,j,k+1) - 30*hxz(i,j,k)   &
                     + 16*hxz(i,j,k-1) -    hxz(i,j,k-2) ) / dzsq12
      d2_hh22(3,3) = (   -hyy(i,j,k+2) + 16*hyy(i,j,k+1) - 30*hyy(i,j,k)   &
                     + 16*hyy(i,j,k-1) -    hyy(i,j,k-2) ) / dzsq12
      d2_hh23(3,3) = (   -hyz(i,j,k+2) + 16*hyz(i,j,k+1) - 30*hyz(i,j,k)   &
                     + 16*hyz(i,j,k-1) -    hyz(i,j,k-2) ) / dzsq12
      d2_hh33(3,3) = (   -hzz(i,j,k+2) + 16*hzz(i,j,k+1) - 30*hzz(i,j,k)   &
                     + 16*hzz(i,j,k-1) -    hzz(i,j,k-2) ) / dzsq12

      d2_hh11(1,2) = (   -hxx(i-2,j+2,k) +  8*hxx(i-1,j+2,k) -  8*hxx(i+1,j+2,k) +   hxx(i+2,j+2,k)   &
                      + 8*hxx(i-2,j+1,k) - 64*hxx(i-1,j+1,k) + 64*hxx(i+1,j+1,k) - 8*hxx(i+2,j+1,k)   &
                      - 8*hxx(i-2,j-1,k) + 64*hxx(i-1,j-1,k) - 64*hxx(i+1,j-1,k) + 8*hxx(i+2,j-1,k)   &
                      +   hxx(i-2,j-2,k) -  8*hxx(i-1,j-2,k) +  8*hxx(i+1,j-2,k) -   hxx(i+2,j-2,k) ) / dxdy144
      d2_hh12(1,2) = (   -hxy(i-2,j+2,k) +  8*hxy(i-1,j+2,k) -  8*hxy(i+1,j+2,k) +   hxy(i+2,j+2,k)   &
                      + 8*hxy(i-2,j+1,k) - 64*hxy(i-1,j+1,k) + 64*hxy(i+1,j+1,k) - 8*hxy(i+2,j+1,k)   &
                      - 8*hxy(i-2,j-1,k) + 64*hxy(i-1,j-1,k) - 64*hxy(i+1,j-1,k) + 8*hxy(i+2,j-1,k)   &
                      +   hxy(i-2,j-2,k) -  8*hxy(i-1,j-2,k) +  8*hxy(i+1,j-2,k) -   hxy(i+2,j-2,k) ) / dxdy144
      d2_hh13(1,2) = (   -hxz(i-2,j+2,k) +  8*hxz(i-1,j+2,k) -  8*hxz(i+1,j+2,k) +   hxz(i+2,j+2,k)   &
                      + 8*hxz(i-2,j+1,k) - 64*hxz(i-1,j+1,k) + 64*hxz(i+1,j+1,k) - 8*hxz(i+2,j+1,k)   &
                      - 8*hxz(i-2,j-1,k) + 64*hxz(i-1,j-1,k) - 64*hxz(i+1,j-1,k) + 8*hxz(i+2,j-1,k)   &
                      +   hxz(i-2,j-2,k) -  8*hxz(i-1,j-2,k) +  8*hxz(i+1,j-2,k) -   hxz(i+2,j-2,k) ) / dxdy144
      d2_hh22(1,2) = (   -hyy(i-2,j+2,k) +  8*hyy(i-1,j+2,k) -  8*hyy(i+1,j+2,k) +   hyy(i+2,j+2,k)   &
                      + 8*hyy(i-2,j+1,k) - 64*hyy(i-1,j+1,k) + 64*hyy(i+1,j+1,k) - 8*hyy(i+2,j+1,k)   &
                      - 8*hyy(i-2,j-1,k) + 64*hyy(i-1,j-1,k) - 64*hyy(i+1,j-1,k) + 8*hyy(i+2,j-1,k)   &
                      +   hyy(i-2,j-2,k) -  8*hyy(i-1,j-2,k) +  8*hyy(i+1,j-2,k) -   hyy(i+2,j-2,k) ) / dxdy144
      d2_hh23(1,2) = (   -hyz(i-2,j+2,k) +  8*hyz(i-1,j+2,k) -  8*hyz(i+1,j+2,k) +   hyz(i+2,j+2,k)   &
                      + 8*hyz(i-2,j+1,k) - 64*hyz(i-1,j+1,k) + 64*hyz(i+1,j+1,k) - 8*hyz(i+2,j+1,k)   &
                      - 8*hyz(i-2,j-1,k) + 64*hyz(i-1,j-1,k) - 64*hyz(i+1,j-1,k) + 8*hyz(i+2,j-1,k)   &
                      +   hyz(i-2,j-2,k) -  8*hyz(i-1,j-2,k) +  8*hyz(i+1,j-2,k) -   hyz(i+2,j-2,k) ) / dxdy144
      d2_hh33(1,2) = (   -hzz(i-2,j+2,k) +  8*hzz(i-1,j+2,k) -  8*hzz(i+1,j+2,k) +   hzz(i+2,j+2,k)   &
                      + 8*hzz(i-2,j+1,k) - 64*hzz(i-1,j+1,k) + 64*hzz(i+1,j+1,k) - 8*hzz(i+2,j+1,k)   &
                      - 8*hzz(i-2,j-1,k) + 64*hzz(i-1,j-1,k) - 64*hzz(i+1,j-1,k) + 8*hzz(i+2,j-1,k)   &
                      +   hzz(i-2,j-2,k) -  8*hzz(i-1,j-2,k) +  8*hzz(i+1,j-2,k) -   hzz(i+2,j-2,k) ) / dxdy144

      d2_hh11(1,3) = (   -hxx(i-2,j,k+2) +  8*hxx(i-1,j,k+2) -  8*hxx(i+1,j,k+2) +   hxx(i+2,j,k+2)   &
                      + 8*hxx(i-2,j,k+1) - 64*hxx(i-1,j,k+1) + 64*hxx(i+1,j,k+1) - 8*hxx(i+2,j,k+1)   &
                      - 8*hxx(i-2,j,k-1) + 64*hxx(i-1,j,k-1) - 64*hxx(i+1,j,k-1) + 8*hxx(i+2,j,k-1)   &
                      +   hxx(i-2,j,k-2) -  8*hxx(i-1,j,k-2) +  8*hxx(i+1,j,k-2) -   hxx(i+2,j,k-2) ) / dxdz144
      d2_hh12(1,3) = (   -hxy(i-2,j,k+2) +  8*hxy(i-1,j,k+2) -  8*hxy(i+1,j,k+2) +   hxy(i+2,j,k+2)   &
                      + 8*hxy(i-2,j,k+1) - 64*hxy(i-1,j,k+1) + 64*hxy(i+1,j,k+1) - 8*hxy(i+2,j,k+1)   &
                      - 8*hxy(i-2,j,k-1) + 64*hxy(i-1,j,k-1) - 64*hxy(i+1,j,k-1) + 8*hxy(i+2,j,k-1)   &
                      +   hxy(i-2,j,k-2) -  8*hxy(i-1,j,k-2) +  8*hxy(i+1,j,k-2) -   hxy(i+2,j,k-2) ) / dxdz144
      d2_hh13(1,3) = (   -hxz(i-2,j,k+2) +  8*hxz(i-1,j,k+2) -  8*hxz(i+1,j,k+2) +   hxz(i+2,j,k+2)   &
                      + 8*hxz(i-2,j,k+1) - 64*hxz(i-1,j,k+1) + 64*hxz(i+1,j,k+1) - 8*hxz(i+2,j,k+1)   &
                      - 8*hxz(i-2,j,k-1) + 64*hxz(i-1,j,k-1) - 64*hxz(i+1,j,k-1) + 8*hxz(i+2,j,k-1)   &
                      +   hxz(i-2,j,k-2) -  8*hxz(i-1,j,k-2) +  8*hxz(i+1,j,k-2) -   hxz(i+2,j,k-2) ) / dxdz144
      d2_hh22(1,3) = (   -hyy(i-2,j,k+2) +  8*hyy(i-1,j,k+2) -  8*hyy(i+1,j,k+2) +   hyy(i+2,j,k+2)   &
                      + 8*hyy(i-2,j,k+1) - 64*hyy(i-1,j,k+1) + 64*hyy(i+1,j,k+1) - 8*hyy(i+2,j,k+1)   &
                      - 8*hyy(i-2,j,k-1) + 64*hyy(i-1,j,k-1) - 64*hyy(i+1,j,k-1) + 8*hyy(i+2,j,k-1)   &
                      +   hyy(i-2,j,k-2) -  8*hyy(i-1,j,k-2) +  8*hyy(i+1,j,k-2) -   hyy(i+2,j,k-2) ) / dxdz144
      d2_hh23(1,3) = (   -hyz(i-2,j,k+2) +  8*hyz(i-1,j,k+2) -  8*hyz(i+1,j,k+2) +   hyz(i+2,j,k+2)   &
                      + 8*hyz(i-2,j,k+1) - 64*hyz(i-1,j,k+1) + 64*hyz(i+1,j,k+1) - 8*hyz(i+2,j,k+1)   &
                      - 8*hyz(i-2,j,k-1) + 64*hyz(i-1,j,k-1) - 64*hyz(i+1,j,k-1) + 8*hyz(i+2,j,k-1)   &
                      +   hyz(i-2,j,k-2) -  8*hyz(i-1,j,k-2) +  8*hyz(i+1,j,k-2) -   hyz(i+2,j,k-2) ) / dxdz144
      d2_hh33(1,3) = (   -hzz(i-2,j,k+2) +  8*hzz(i-1,j,k+2) -  8*hzz(i+1,j,k+2) +   hzz(i+2,j,k+2)   &
                      + 8*hzz(i-2,j,k+1) - 64*hzz(i-1,j,k+1) + 64*hzz(i+1,j,k+1) - 8*hzz(i+2,j,k+1)   &
                      - 8*hzz(i-2,j,k-1) + 64*hzz(i-1,j,k-1) - 64*hzz(i+1,j,k-1) + 8*hzz(i+2,j,k-1)   &
                      +   hzz(i-2,j,k-2) -  8*hzz(i-1,j,k-2) +  8*hzz(i+1,j,k-2) -   hzz(i+2,j,k-2) ) / dxdz144

      d2_hh11(2,3) = (   -hxx(i,j-2,k+2) +  8*hxx(i,j-1,k+2) -  8*hxx(i,j+1,k+2) +   hxx(i,j+2,k+2)   &
                      + 8*hxx(i,j-2,k+1) - 64*hxx(i,j-1,k+1) + 64*hxx(i,j+1,k+1) - 8*hxx(i,j+2,k+1)   &
                      - 8*hxx(i,j-2,k-1) + 64*hxx(i,j-1,k-1) - 64*hxx(i,j+1,k-1) + 8*hxx(i,j+2,k-1)   &
                      +   hxx(i,j-2,k-2) -  8*hxx(i,j-1,k-2) +  8*hxx(i,j+1,k-2) -   hxx(i,j+2,k-2) ) / dydz144
      d2_hh12(2,3) = (   -hxy(i,j-2,k+2) +  8*hxy(i,j-1,k+2) -  8*hxy(i,j+1,k+2) +   hxy(i,j+2,k+2)   &
                      + 8*hxy(i,j-2,k+1) - 64*hxy(i,j-1,k+1) + 64*hxy(i,j+1,k+1) - 8*hxy(i,j+2,k+1)   &
                      - 8*hxy(i,j-2,k-1) + 64*hxy(i,j-1,k-1) - 64*hxy(i,j+1,k-1) + 8*hxy(i,j+2,k-1)   &
                      +   hxy(i,j-2,k-2) -  8*hxy(i,j-1,k-2) +  8*hxy(i,j+1,k-2) -   hxy(i,j+2,k-2) ) / dydz144
      d2_hh13(2,3) = (   -hxz(i,j-2,k+2) +  8*hxz(i,j-1,k+2) -  8*hxz(i,j+1,k+2) +   hxz(i,j+2,k+2)   &
                      + 8*hxz(i,j-2,k+1) - 64*hxz(i,j-1,k+1) + 64*hxz(i,j+1,k+1) - 8*hxz(i,j+2,k+1)   &
                      - 8*hxz(i,j-2,k-1) + 64*hxz(i,j-1,k-1) - 64*hxz(i,j+1,k-1) + 8*hxz(i,j+2,k-1)   &
                      +   hxz(i,j-2,k-2) -  8*hxz(i,j-1,k-2) +  8*hxz(i,j+1,k-2) -   hxz(i,j+2,k-2) ) / dydz144
      d2_hh22(2,3) = (   -hyy(i,j-2,k+2) +  8*hyy(i,j-1,k+2) -  8*hyy(i,j+1,k+2) +   hyy(i,j+2,k+2)   &
                      + 8*hyy(i,j-2,k+1) - 64*hyy(i,j-1,k+1) + 64*hyy(i,j+1,k+1) - 8*hyy(i,j+2,k+1)   &
                      - 8*hyy(i,j-2,k-1) + 64*hyy(i,j-1,k-1) - 64*hyy(i,j+1,k-1) + 8*hyy(i,j+2,k-1)   &
                      +   hyy(i,j-2,k-2) -  8*hyy(i,j-1,k-2) +  8*hyy(i,j+1,k-2) -   hyy(i,j+2,k-2) ) / dydz144
      d2_hh23(2,3) = (   -hyz(i,j-2,k+2) +  8*hyz(i,j-1,k+2) -  8*hyz(i,j+1,k+2) +   hyz(i,j+2,k+2)   &
                      + 8*hyz(i,j-2,k+1) - 64*hyz(i,j-1,k+1) + 64*hyz(i,j+1,k+1) - 8*hyz(i,j+2,k+1)   &
                      - 8*hyz(i,j-2,k-1) + 64*hyz(i,j-1,k-1) - 64*hyz(i,j+1,k-1) + 8*hyz(i,j+2,k-1)   &
                      +   hyz(i,j-2,k-2) -  8*hyz(i,j-1,k-2) +  8*hyz(i,j+1,k-2) -   hyz(i,j+2,k-2) ) / dydz144
      d2_hh33(2,3) = (   -hzz(i,j-2,k+2) +  8*hzz(i,j-1,k+2) -  8*hzz(i,j+1,k+2) +   hzz(i,j+2,k+2)   &
                      + 8*hzz(i,j-2,k+1) - 64*hzz(i,j-1,k+1) + 64*hzz(i,j+1,k+1) - 8*hzz(i,j+2,k+1)   &
                      - 8*hzz(i,j-2,k-1) + 64*hzz(i,j-1,k-1) - 64*hzz(i,j+1,k-1) + 8*hzz(i,j+2,k-1)   &
                      +   hzz(i,j-2,k-2) -  8*hzz(i,j-1,k-2) +  8*hzz(i,j+1,k-2) -   hzz(i,j+2,k-2) ) / dydz144

      d2_hh11(2,1) = d2_hh11(1,2)
      d2_hh12(2,1) = d2_hh12(1,2)
      d2_hh13(2,1) = d2_hh13(1,2)
      d2_hh22(2,1) = d2_hh22(1,2)
      d2_hh23(2,1) = d2_hh23(1,2)
      d2_hh33(2,1) = d2_hh33(1,2)

      d2_hh11(3,1) = d2_hh11(1,3)
      d2_hh12(3,1) = d2_hh12(1,3)
      d2_hh13(3,1) = d2_hh13(1,3)
      d2_hh22(3,1) = d2_hh22(1,3)
      d2_hh23(3,1) = d2_hh23(1,3)
      d2_hh33(3,1) = d2_hh33(1,3)

      d2_hh11(3,2) = d2_hh11(2,3)
      d2_hh12(3,2) = d2_hh12(2,3)
      d2_hh13(3,2) = d2_hh13(2,3)
      d2_hh22(3,2) = d2_hh22(2,3)
      d2_hh23(3,2) = d2_hh23(2,3)
      d2_hh33(3,2) = d2_hh33(2,3)

    end if
    

    if (use_jacobian) then
      call LeanBSSN_apply_jacobian(d1_trk, jac)
      call LeanBSSN_apply_jacobian(d1_gammat1, jac)
      call LeanBSSN_apply_jacobian(d1_gammat2, jac)
      call LeanBSSN_apply_jacobian(d1_gammat3, jac)
      call LeanBSSN_apply_jacobian(d1_aa11, jac)
      call LeanBSSN_apply_jacobian(d1_aa12, jac)
      call LeanBSSN_apply_jacobian(d1_aa13, jac)
      call LeanBSSN_apply_jacobian(d1_aa22, jac)
      call LeanBSSN_apply_jacobian(d1_aa23, jac)
      call LeanBSSN_apply_jacobian(d1_aa33, jac)

      call LeanBSSN_apply_jacobian2(d1_ww, d2_ww, jac, hes)
      call LeanBSSN_apply_jacobian2(d1_hh11, d2_hh11, jac, hes)
      call LeanBSSN_apply_jacobian2(d1_hh12, d2_hh12, jac, hes)
      call LeanBSSN_apply_jacobian2(d1_hh13, d2_hh13, jac, hes)
      call LeanBSSN_apply_jacobian2(d1_hh22, d2_hh22, jac, hes)
      call LeanBSSN_apply_jacobian2(d1_hh23, d2_hh23, jac, hes)
      call LeanBSSN_apply_jacobian2(d1_hh33, d2_hh33, jac, hes)
    end if

    d1_gammat(1,:) = d1_gammat1(:)
    d1_gammat(2,:) = d1_gammat2(:)
    d1_gammat(3,:) = d1_gammat3(:)


    d1_hh(1,1,:) = d1_hh11(:)
    d1_hh(1,2,:) = d1_hh12(:)
    d1_hh(1,3,:) = d1_hh13(:)
    d1_hh(2,2,:) = d1_hh22(:)
    d1_hh(2,3,:) = d1_hh23(:)
    d1_hh(3,3,:) = d1_hh33(:)
    d1_hh(2,1,:) = d1_hh(1,2,:)
    d1_hh(3,1,:) = d1_hh(1,3,:)
    d1_hh(3,2,:) = d1_hh(2,3,:)

    d1_aa(1,1,:) = d1_aa11(:)
    d1_aa(1,2,:) = d1_aa12(:)
    d1_aa(1,3,:) = d1_aa13(:)
    d1_aa(2,2,:) = d1_aa22(:)
    d1_aa(2,3,:) = d1_aa23(:)
    d1_aa(3,3,:) = d1_aa33(:)
    d1_aa(2,1,:) = d1_aa(1,2,:)
    d1_aa(3,1,:) = d1_aa(1,3,:)
    d1_aa(3,2,:) = d1_aa(2,3,:)

    d2_hh(1,1,:,:) = d2_hh11(:,:)
    d2_hh(1,2,:,:) = d2_hh12(:,:)
    d2_hh(1,3,:,:) = d2_hh13(:,:)
    d2_hh(2,2,:,:) = d2_hh22(:,:)
    d2_hh(2,3,:,:) = d2_hh23(:,:)
    d2_hh(3,3,:,:) = d2_hh33(:,:)
    d2_hh(2,1,:,:) = d2_hh(1,2,:,:)
    d2_hh(3,1,:,:) = d2_hh(1,3,:,:)
    d2_hh(3,2,:,:) = d2_hh(2,3,:,:)
    
    cf1 = 0
    do a = 1, 3
      do b = 1, 3
        do c = b, 3
          cf1(a,b,c) = 0.5d0 * (d1_hh(a,b,c) + d1_hh(a,c,b) - d1_hh(b,c,a))
        end do
      end do
    end do
    cf1(:,2,1) = cf1(:,1,2)
    cf1(:,3,1) = cf1(:,1,3)
    cf1(:,3,2) = cf1(:,2,3)

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


    
    cd2_ww = d2_ww
    cd1_aa   = d1_aa
    do a = 1, 3
      do b = a, 3
        do l = 1, 3
          cd2_ww(a,b) = cd2_ww(a,b) - cf2(l,a,b) * d1_ww(l)
        do m = 1, 3
            cd1_aa(a,b,l) = cd1_aa(a,b,l) - cf2(m,a,l) * aa(b,m) - cf2(m,b,l) * aa(a,m)
        end do
          end do
        end do
      end do
    cd2_ww(2,1)   = cd2_ww(1,2)
    cd2_ww(3,1)   = cd2_ww(1,3)
    cd2_ww(3,2)   = cd2_ww(2,3)
    cd1_aa(2,1,:) = cd1_aa(1,2,:)
    cd1_aa(3,1,:) = cd1_aa(1,3,:)
    cd1_aa(3,2,:) = cd1_aa(2,3,:)
    


    
    
    ri_1 = 0
    ri_2 = 0
    ri_3 = 0
    c_ri_ww = 0
    c_ri_hh = 0

    tr_cd2_ww = 0
    tr_dww_dww = 0
    do l = 1, 3
      do m = 1, 3
        tr_cd2_ww  = tr_cd2_ww  + hu(l,m) * cd2_ww(l,m)
        tr_dww_dww = tr_dww_dww + hu(l,m) * d1_ww(l) * d1_ww(m)
      end do
    end do
    
    c_ri_ww = ww * ( cd2_ww + hh * tr_cd2_ww ) - 2 * hh * tr_dww_dww


    do a = 1, 3
      do b = a, 3
        do l = 1, 3
          ri_1(a,b) = ri_1(a,b) + hh(l,a) * d1_gammat(l,b) / 2              &
                                + hh(l,b) * d1_gammat(l,a) / 2              &
                                + gammat(l) * (cf1(a,b,l) + cf1(b,a,l)) / 2
          do m = 1, 3
            ri_2(a,b) = ri_2(a,b) - hu(l,m) * d2_hh(a,b,l,m) / 2
            do n = 1, 3
              ri_3(a,b) = ri_3(a,b) + hu(l,m)                              &
                                    * ( cf2(n,l,a) * cf1(b,n,m)            &
                                    + cf2(n,l,b) * (cf1(a,n,m) + cf1(n,m,a)) )
            end do
          end do
        end do
      end do
    end do
    
    c_ri_hh = ww*ww * (ri_1 + ri_2 + ri_3)
    c_ri    = c_ri_ww + c_ri_hh

    c_ri(2,1) = c_ri(1,2)
    c_ri(3,1) = c_ri(1,3)
    c_ri(3,2) = c_ri(2,3)

    trr = 0
    do m = 1, 3
      do n = 1, 3
        trr = trr + hu(m,n) * c_ri(m,n)
      end do
    end do
    


    
    sq_aa = 0
    a2    = 0
    do m = 1, 3
      do n = 1, 3
        do p = 1, 3
          do q = 1, 3
            a2(m,n) = a2(m,n) + hu(p,q) * aa(m,p) * aa(n,q)
          end do
        end do
        sq_aa = sq_aa + hu(m,n) * a2(m,n)
      end do
    end do
    


    
    ham = trr + 2 * trk**2 / 3 - sq_aa

    mom = -2 * d1_trk / 3
    do a = 1, 3
      do l = 1, 3
        do m = 1, 3
          mom(a) = mom(a) + hu(l,m) * ( cd1_aa(a,l,m) - 3 * aa(a,l) * d1_ww(m) / ww )
        end do
      end do
    end do
    


    
    
    
    
    
    
    
    
    
    
    Tab = 0
    if (stress_energy_state /= 0) then
       Tab(4,4) = eTtt(i,j,k)
       Tab(4,1) = eTtx(i,j,k)
       Tab(4,2) = eTty(i,j,k)
       Tab(4,3) = eTtz(i,j,k)
       Tab(1,1) = eTxx(i,j,k)
       Tab(1,2) = eTxy(i,j,k)
       Tab(1,3) = eTxz(i,j,k)
       Tab(2,2) = eTyy(i,j,k)
       Tab(2,3) = eTyz(i,j,k)
       Tab(3,3) = eTzz(i,j,k)
       Tab(1,4) = Tab(4,1)
       Tab(2,4) = Tab(4,2)
       Tab(3,4) = Tab(4,3)
       Tab(2,1) = Tab(1,2)
       Tab(3,1) = Tab(1,3)
       Tab(3,2) = Tab(2,3)

       srcE = Tab(4,4)
       do m = 1, 3
          srcE = srcE - 2 * beta(m) * Tab(m,4)
          do n = 1, 3
             srcE = srcE + beta(m) * beta(n) * Tab(m,n)
          end do
       end do
       srcE = srcE / (alph * alph)

       srcjdi = 0
       do a = 1, 3
          do m = 1, 3
             srcjdi(a) = srcjdi(a) + beta(m) * Tab(a,m)
          end do
       end do
       srcjdi = (srcjdi - Tab(1:3,4)) / alph

       
       ham = ham - pi16 * srcE
       mom = mom - pi8  * srcjdi

    end if

    
    hc(i,j,k)  = ham

    mcx(i,j,k) = mom(1)
    mcy(i,j,k) = mom(2)
    mcz(i,j,k) = mom(3)
    

  end do
  end do
  end do
  !$OMP END PARALLEL DO

end subroutine LeanBSSN_bssn_constraints
