





















































































































 


































































































































































 



































































































































































































 



















































 










 








































































































































































































































 










                                
                                

                                
                                

                                
                                

                                




















































































subroutine Proca_calc_rhs( cctk_dim,cctk_gsh,cctk_lsh,cctk_lbnd,cctk_ubnd,cctk_tile_min,cctk_tile_max,cctk_ash,cctk_alignment,cctk_&
  &alignment_offset,cctk_from,cctk_to,cctk_bbox,cctk_delta_time,cctk_time,cctk_delta_space,cctk_origin_space,cctk_levfac,cctk_levof&
  &f,cctk_levoffdenom,cctk_timefac,cctk_convlevel,cctk_convfac,cctk_nghostzones,cctk_iteration,cctkGH,cctk_ash1,cctk_ash2,cctk_ash3&
  &, Aphi,Aphi_p,Aphi_p_p,Ax,Ax_p,Ax_p_p,Ay,Ay_p,Ay_p_p,Az,Az_p,Az_p_p,Ex,Ex_p,Ex_p_p,Ey,Ey_p,Ey_p_p,Ez,Ez_p,Ez_p_p,Zeta,Zeta_p,Zet&
  &a_p_p,alp,alp_p,alp_p_p,betax,betax_p,betax_p_p,betay,betay_p,betay_p_p,betaz,betaz_p,betaz_p_p,chi,coarse_dx,coarse_dy,coarse_d&
  &z,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,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,hxx,hxy,hxz,hyy,hyz,hzz,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,psi,psix,psixx,psixy,psixz,psi&
  &y,psiyy,psiyz,psiz,psizz,r,rhs_Aphi,rhs_Ax,rhs_Ay,rhs_Az,rhs_Ex,rhs_Ey,rhs_Ez,rhs_Zeta,shift_state,stress_energy_state,tracek,x,&
  &y,z )

  implicit none

  INTEGER cctk_dim
  integer, parameter :: cctki_use_cctk_dim = kind(cctk_dim)
  INTEGER cctk_gsh (cctk_dim)
  integer, parameter :: cctki_use_cctk_gsh = kind(cctk_gsh)
  INTEGER cctk_lsh (cctk_dim)
  integer, parameter :: cctki_use_cctk_lsh = kind(cctk_lsh)
  INTEGER cctk_lbnd (cctk_dim)
  integer, parameter :: cctki_use_cctk_lbnd = kind(cctk_lbnd)
  INTEGER cctk_ubnd (cctk_dim)
  integer, parameter :: cctki_use_cctk_ubnd = kind(cctk_ubnd)
  INTEGER cctk_tile_min (cctk_dim)
  integer, parameter :: cctki_use_cctk_tile_min = kind(cctk_tile_min)
  INTEGER cctk_tile_max (cctk_dim)
  integer, parameter :: cctki_use_cctk_tile_max = kind(cctk_tile_max)
  INTEGER cctk_ash (cctk_dim)
  integer, parameter :: cctki_use_cctk_ash = kind(cctk_ash)
  INTEGER cctk_alignment
  integer, parameter :: cctki_use_cctk_alignment = kind(cctk_alignment)
  INTEGER cctk_alignment_offset
  integer, parameter :: cctki_use_cctk_alignment_offset = kind(cctk_alignment_offset)
  INTEGER cctk_from (cctk_dim)
  integer, parameter :: cctki_use_cctk_from = kind(cctk_from)
  INTEGER cctk_to (cctk_dim)
  integer, parameter :: cctki_use_cctk_to = kind(cctk_to)
  INTEGER cctk_bbox (2*cctk_dim)
  integer, parameter :: cctki_use_cctk_bbox = kind(cctk_bbox)
  REAL*8 cctk_delta_time
  integer, parameter :: cctki_use_cctk_delta_time = kind(cctk_delta_time)
  REAL*8 cctk_time
  integer, parameter :: cctki_use_cctk_time = kind(cctk_time)
  REAL*8 cctk_delta_space (cctk_dim)
  integer, parameter :: cctki_use_cctk_delta_space = kind(cctk_delta_space)
  REAL*8 cctk_origin_space (cctk_dim)
  integer, parameter :: cctki_use_cctk_origin_space = kind(cctk_origin_space)
  INTEGER cctk_levfac (cctk_dim)
  integer, parameter :: cctki_use_cctk_levfac = kind(cctk_levfac)
  INTEGER cctk_levoff (cctk_dim)
  integer, parameter :: cctki_use_cctk_levoff = kind(cctk_levoff)
  INTEGER cctk_levoffdenom (cctk_dim)
  integer, parameter :: cctki_use_cctk_levoffdenom = kind(cctk_levoffdenom)
  INTEGER cctk_timefac
  integer, parameter :: cctki_use_cctk_timefac = kind(cctk_timefac)
  INTEGER cctk_convlevel
  integer, parameter :: cctki_use_cctk_convlevel = kind(cctk_convlevel)
  INTEGER cctk_convfac
  integer, parameter :: cctki_use_cctk_convfac = kind(cctk_convfac)
  INTEGER cctk_nghostzones (cctk_dim)
  integer, parameter :: cctki_use_cctk_nghostzones = kind(cctk_nghostzones)
  INTEGER cctk_iteration
  integer, parameter :: cctki_use_cctk_iteration = kind(cctk_iteration)
  integer*8 cctkGH
  integer, parameter :: cctki_use_cctkGH = kind(cctkGH)
  INTEGER cctk_ash1
  integer, parameter :: cctki_use_cctk_ash1 = kind(cctk_ash1)
  INTEGER cctk_ash2
  integer, parameter :: cctki_use_cctk_ash2 = kind(cctk_ash2)
  INTEGER cctk_ash3
  integer, parameter :: cctki_use_cctk_ash3 = kind(cctk_ash3)
  REAL*8 Aphi (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Aphi = kind(Aphi)
  REAL*8 Aphi_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Aphi_p = kind(Aphi_p)
  REAL*8 Aphi_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Aphi_p_p = kind(Aphi_p_p)
  REAL*8 Ax (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Ax = kind(Ax)
  REAL*8 Ax_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Ax_p = kind(Ax_p)
  REAL*8 Ax_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Ax_p_p = kind(Ax_p_p)
  REAL*8 Ay (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Ay = kind(Ay)
  REAL*8 Ay_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Ay_p = kind(Ay_p)
  REAL*8 Ay_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Ay_p_p = kind(Ay_p_p)
  REAL*8 Az (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Az = kind(Az)
  REAL*8 Az_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Az_p = kind(Az_p)
  REAL*8 Az_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Az_p_p = kind(Az_p_p)
  REAL*8 Ex (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Ex = kind(Ex)
  REAL*8 Ex_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Ex_p = kind(Ex_p)
  REAL*8 Ex_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Ex_p_p = kind(Ex_p_p)
  REAL*8 Ey (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Ey = kind(Ey)
  REAL*8 Ey_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Ey_p = kind(Ey_p)
  REAL*8 Ey_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Ey_p_p = kind(Ey_p_p)
  REAL*8 Ez (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Ez = kind(Ez)
  REAL*8 Ez_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Ez_p = kind(Ez_p)
  REAL*8 Ez_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Ez_p_p = kind(Ez_p_p)
  REAL*8 Zeta (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Zeta = kind(Zeta)
  REAL*8 Zeta_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Zeta_p = kind(Zeta_p)
  REAL*8 Zeta_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Zeta_p_p = kind(Zeta_p_p)
  REAL*8 alp (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_alp = kind(alp)
  REAL*8 alp_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_alp_p = kind(alp_p)
  REAL*8 alp_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_alp_p_p = kind(alp_p_p)
  REAL*8 betax (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_betax = kind(betax)
  REAL*8 betax_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_betax_p = kind(betax_p)
  REAL*8 betax_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_betax_p_p = kind(betax_p_p)
  REAL*8 betay (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_betay = kind(betay)
  REAL*8 betay_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_betay_p = kind(betay_p)
  REAL*8 betay_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_betay_p_p = kind(betay_p_p)
  REAL*8 betaz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_betaz = kind(betaz)
  REAL*8 betaz_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_betaz_p = kind(betaz_p)
  REAL*8 betaz_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_betaz_p_p = kind(betaz_p_p)
  REAL*8 chi (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_chi = kind(chi)
  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)
  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 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 hxx (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_hxx = kind(hxx)
  REAL*8 hxy (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_hxy = kind(hxy)
  REAL*8 hxz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_hxz = kind(hxz)
  REAL*8 hyy (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_hyy = kind(hyy)
  REAL*8 hyz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_hyz = kind(hyz)
  REAL*8 hzz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_hzz = kind(hzz)
  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 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_Aphi (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_Aphi = kind(rhs_Aphi)
  REAL*8 rhs_Ax (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_Ax = kind(rhs_Ax)
  REAL*8 rhs_Ay (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_Ay = kind(rhs_Ay)
  REAL*8 rhs_Az (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_Az = kind(rhs_Az)
  REAL*8 rhs_Ex (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_Ex = kind(rhs_Ex)
  REAL*8 rhs_Ey (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_Ey = kind(rhs_Ey)
  REAL*8 rhs_Ez (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_Ez = kind(rhs_Ez)
  REAL*8 rhs_Zeta (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_Zeta = kind(rhs_Zeta)
  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 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 ProcaEvolve_MaxNumConstrainedVars
  integer, parameter :: cctki_use_ProcaEvolve_MaxNumConstrainedVars = kind(ProcaEvolve_MaxNumConstrainedVars)
  INTEGER*4 ProcaEvolve_MaxNumEvolvedVars
  integer, parameter :: cctki_use_ProcaEvolve_MaxNumEvolvedVars = kind(ProcaEvolve_MaxNumEvolvedVars)
  INTEGER*4 ProcaEvolve_MaxNumSandRVars
  integer, parameter :: cctki_use_ProcaEvolve_MaxNumSandRVars = kind(ProcaEvolve_MaxNumSandRVars)
  COMMON /ProcaEvolverest/ProcaEvolve_MaxNumConstrainedVars, ProcaEvolve_MaxNumEvolvedVars, ProcaEvolve_MaxNumSandRVars
  REAL*8 Zeta_Omega_fac
  integer, parameter :: cctki_use_Zeta_Omega_fac = kind(Zeta_Omega_fac)
  REAL*8 chi_floor
  integer, parameter :: cctki_use_chi_floor = kind(chi_floor)
  REAL*8 conf_fac_exponent
  integer, parameter :: cctki_use_conf_fac_exponent = kind(conf_fac_exponent)
  REAL*8 kappa
  integer, parameter :: cctki_use_kappa = kind(kappa)
  INTEGER*4 derivs_order
  integer, parameter :: cctki_use_derivs_order = kind(derivs_order)
  INTEGER*4 use_advection_stencils
  integer, parameter :: cctki_use_use_advection_stencils = kind(use_advection_stencils)
  COMMON /ProcaEvolvepriv/Zeta_Omega_fac, chi_floor, conf_fac_exponent, kappa, derivs_order, use_advection_stencils
  INTEGER*4 CCTKH0
  integer, parameter :: cctki_use_CCTKH0 = kind(CCTKH0)
  INTEGER*4 CCTKH1
  integer, parameter :: cctki_use_CCTKH1 = kind(CCTKH1)
  INTEGER*4 CCTKH2
  integer, parameter :: cctki_use_CCTKH2 = kind(CCTKH2)
  INTEGER*4 CCTKH3
  integer, parameter :: cctki_use_CCTKH3 = kind(CCTKH3)
  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 CCTKH4
  integer, parameter :: cctki_use_CCTKH4 = kind(CCTKH4)
  INTEGER*4 MoL_Num_SaveAndRestore_Vars
  integer, parameter :: cctki_use_MoL_Num_SaveAndRestore_Vars = kind(MoL_Num_SaveAndRestore_Vars)
  INTEGER*4 CCTKH5
  integer, parameter :: cctki_use_CCTKH5 = kind(CCTKH5)
  COMMON /METHODOFLINESrest/CCTKH0, CCTKH1, CCTKH2, CCTKH3, MoL_Num_Constrained_Vars, MoL_Num_Evolved_Vars, CCTKH4, MoL_Num_SaveAnd&
  &Restore_Vars, CCTKH5
  REAL*8 A0 (3)
  integer, parameter :: cctki_use_A0 = kind(A0)
  REAL*8 Aphi0
  integer, parameter :: cctki_use_Aphi0 = kind(Aphi0)
  REAL*8 E0 (3)
  integer, parameter :: cctki_use_E0 = kind(E0)
  REAL*8 mu
  integer, parameter :: cctki_use_mu = kind(mu)
  integer*8 evolution_method
  integer, parameter :: cctki_use_evolution_method = kind(evolution_method)
  integer*8 CCTKH6
  integer, parameter :: cctki_use_CCTKH6 = kind(CCTKH6)
  INTEGER*4 n_A (3)
  integer, parameter :: cctki_use_n_A = kind(n_A)
  INTEGER*4 n_Aphi
  integer, parameter :: cctki_use_n_Aphi = kind(n_Aphi)
  INTEGER*4 n_E (3)
  integer, parameter :: cctki_use_n_E = kind(n_E)
  INTEGER*4 n_Zeta
  integer, parameter :: cctki_use_n_Zeta = kind(n_Zeta)
  COMMON /PROCABASErest/A0, Aphi0, E0, mu, evolution_method, CCTKH6, n_A, n_Aphi, n_E, n_Zeta

  
  REAL*8                alph, beta(3)
  REAL*8                hh(3,3), hu(3,3), trk, dethh, ch
  REAL*8                lE(3), lA(3), lAphi, lZeta

  
  REAL*8                d1_alph(3), d1_beta(3,3)
  REAL*8                d1_hh(3,3,3), d1_ch(3)
  REAL*8                d1_lE(3,3), d1_lA(3,3), d1_lZeta(3), d1_lAphi(3)

  
  REAL*8                d2_lA(3,3,3)

  
  REAL*8                ad1_lE(3), ad1_lA(3), ad1_lZeta, ad1_lAphi
  REAL*8                d1_f(3)   

  
  REAL*8                cf1(3,3,3), cf2(3,3,3)

  
  REAL*8                cd_lA(3,3), cd_dA(3,3,3)

  
  REAL*8                rhs_lE(3), rhs_lA(3), rhs_lZeta, rhs_lAphi

  
  REAL*8                dx12, dy12, dz12, dxsq12, dysq12, dzsq12,         &
                           dxdy144, dxdz144, dydz144
  REAL*8                odx60, ody60, odz60, odxsq180, odysq180, odzsq180,&
                           odxdy3600, odxdz3600, odydz3600
  INTEGER*4                 i, j, k
  INTEGER*4                 di, dj, dk
  REAL*8, parameter ::  one = 1
  INTEGER*4                 a, b, c, m, n


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


  
  call Proca_adm2bssn(cctk_dim,cctk_gsh,cctk_lsh,cctk_lbnd,cctk_ubnd,cctk_tile_min,cctk_tile_max,cctk_ash,cctk_alignment,cctk_align&
  &ment_offset,cctk_from,cctk_to,cctk_bbox,cctk_delta_time,cctk_time,cctk_delta_space,cctk_origin_space,cctk_levfac,cctk_levoff,cct&
  &k_levoffdenom,cctk_timefac,cctk_convlevel,cctk_convfac,cctk_nghostzones,cctk_iteration,cctkGH,cctk_ash1,cctk_ash2,cctk_ash3, Aph&
  &i,Aphi_p,Aphi_p_p,Ax,Ax_p,Ax_p_p,Ay,Ay_p,Ay_p_p,Az,Az_p,Az_p_p,Ex,Ex_p,Ex_p_p,Ey,Ey_p,Ey_p_p,Ez,Ez_p,Ez_p_p,Zeta,Zeta_p,Zeta_p_p&
  &,alp,alp_p,alp_p_p,betax,betax_p,betax_p_p,betay,betay_p,betay_p_p,betaz,betaz_p,betaz_p_p,chi,coarse_dx,coarse_dy,coarse_dz,con&
  &formal_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,d&
  &tlapse_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,gxx,gxx_p,gxx_p_p,g&
  &xy,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,hxx,hxy,hxz,hyy,hyz,hzz,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,psi,psix,psixx,psixy,psixz,psiy,psi&
  &yy,psiyz,psiz,psizz,r,rhs_Aphi,rhs_Ax,rhs_Ay,rhs_Az,rhs_Ex,rhs_Ey,rhs_Ez,rhs_Zeta,shift_state,stress_energy_state,tracek,x,y,z)


  !$OMP PARALLEL DO COLLAPSE(3) &
  !$OMP PRIVATE(alph, beta, hh, hu, trk, dethh, ch,&
  !$OMP lE, lA, lAphi, lZeta,&
  !$OMP d1_alph, d1_beta, d1_hh, d1_ch,&
  !$OMP d1_lE, d1_lA, d1_lZeta, d1_lAphi,&
  !$OMP d2_lA, ad1_lE, ad1_lA, ad1_lZeta, ad1_lAphi,&
  !$OMP d1_f, cf1, cf2, cd_lA, cd_dA,&
  !$OMP rhs_lE, rhs_lA, rhs_lZeta, rhs_lAphi,&
  !$OMP i, j, k,&
  !$OMP di, dj, dk,&
  !$OMP a, b, c, m, n)
  do k = 1+cctk_nghostzones(3), cctk_lsh(3)-cctk_nghostzones(3)
  do j = 1+cctk_nghostzones(2), cctk_lsh(2)-cctk_nghostzones(2)
  do i = 1+cctk_nghostzones(1), cctk_lsh(1)-cctk_nghostzones(1)

    
    ch        = chi(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)

    alph      = alp(i,j,k)

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

    lE(1)     = Ex(i,j,k)
    lE(2)     = Ey(i,j,k)
    lE(3)     = Ez(i,j,k)

    lA(1)     = Ax(i,j,k)
    lA(2)     = Ay(i,j,k)
    lA(3)     = Az(i,j,k)

    lZeta     = Zeta(i,j,k)
    lAphi     = Aphi(i,j,k)
    


    
    
    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 == 4) then

      

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

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

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


      
      d1_hh(1,1,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_hh(1,2,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_hh(1,3,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_hh(2,2,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_hh(2,3,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_hh(3,3,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_hh(1,1,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_hh(1,2,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_hh(1,3,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_hh(2,2,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_hh(2,3,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_hh(3,3,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_hh(1,1,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_hh(1,2,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_hh(1,3,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_hh(2,2,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_hh(2,3,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_hh(3,3,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_hh(2,1,:) = d1_hh(1,2,:)
      d1_hh(3,1,:) = d1_hh(1,3,:)
      d1_hh(3,2,:) = d1_hh(2,3,:)


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

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

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

      
      d1_beta(1,1)  = (   -betax(i+2,j,k) + 8*betax(i+1,j,k)                 &
                       - 8*betax(i-1,j,k) +   betax(i-2,j,k) ) / dx12
      d1_beta(2,1)  = (   -betay(i+2,j,k) + 8*betay(i+1,j,k)                 &
                       - 8*betay(i-1,j,k) +   betay(i-2,j,k) ) / dx12
      d1_beta(3,1)  = (   -betaz(i+2,j,k) + 8*betaz(i+1,j,k)                 &
                       - 8*betaz(i-1,j,k) +   betaz(i-2,j,k) ) / dx12

      d1_beta(1,2)  = (   -betax(i,j+2,k) + 8*betax(i,j+1,k)                 &
                       - 8*betax(i,j-1,k) +   betax(i,j-2,k) ) / dy12
      d1_beta(2,2)  = (   -betay(i,j+2,k) + 8*betay(i,j+1,k)                 &
                       - 8*betay(i,j-1,k) +   betay(i,j-2,k) ) / dy12
      d1_beta(3,2)  = (   -betaz(i,j+2,k) + 8*betaz(i,j+1,k)                 &
                       - 8*betaz(i,j-1,k) +   betaz(i,j-2,k) ) / dy12

      d1_beta(1,3)  = (   -betax(i,j,k+2) + 8*betax(i,j,k+1)                 &
                       - 8*betax(i,j,k-1) +   betax(i,j,k-2) ) / dz12
      d1_beta(2,3)  = (   -betay(i,j,k+2) + 8*betay(i,j,k+1)                 &
                       - 8*betay(i,j,k-1) +   betay(i,j,k-2) ) / dz12
      d1_beta(3,3)  = (   -betaz(i,j,k+2) + 8*betaz(i,j,k+1)                 &
                       - 8*betaz(i,j,k-1) +   betaz(i,j,k-2) ) / dz12

      
      d1_lE(1,1) = (   -Ex(i+2,j,k) + 8*Ex(i+1,j,k)               &
                    - 8*Ex(i-1,j,k) +   Ex(i-2,j,k) ) / dx12
      d1_lE(2,1) = (   -Ey(i+2,j,k) + 8*Ey(i+1,j,k)               &
                    - 8*Ey(i-1,j,k) +   Ey(i-2,j,k) ) / dx12
      d1_lE(3,1) = (   -Ez(i+2,j,k) + 8*Ez(i+1,j,k)               &
                    - 8*Ez(i-1,j,k) +   Ez(i-2,j,k) ) / dx12

      d1_lE(1,2) = (   -Ex(i,j+2,k) + 8*Ex(i,j+1,k)               &
                    - 8*Ex(i,j-1,k) +   Ex(i,j-2,k) ) / dy12
      d1_lE(2,2) = (   -Ey(i,j+2,k) + 8*Ey(i,j+1,k)               &
                    - 8*Ey(i,j-1,k) +   Ey(i,j-2,k) ) / dy12
      d1_lE(3,2) = (   -Ez(i,j+2,k) + 8*Ez(i,j+1,k)               &
                    - 8*Ez(i,j-1,k) +   Ez(i,j-2,k) ) / dy12

      d1_lE(1,3) = (   -Ex(i,j,k+2) + 8*Ex(i,j,k+1)               &
                    - 8*Ex(i,j,k-1) +   Ex(i,j,k-2) ) / dz12
      d1_lE(2,3) = (   -Ey(i,j,k+2) + 8*Ey(i,j,k+1)               &
                    - 8*Ey(i,j,k-1) +   Ey(i,j,k-2) ) / dz12
      d1_lE(3,3) = (   -Ez(i,j,k+2) + 8*Ez(i,j,k+1)               &
                    - 8*Ez(i,j,k-1) +   Ez(i,j,k-2) ) / dz12

      
      d1_lA(1,1) = (   -Ax(i+2,j,k) + 8*Ax(i+1,j,k)               &
                    - 8*Ax(i-1,j,k) +   Ax(i-2,j,k) ) / dx12
      d1_lA(2,1) = (   -Ay(i+2,j,k) + 8*Ay(i+1,j,k)               &
                    - 8*Ay(i-1,j,k) +   Ay(i-2,j,k) ) / dx12
      d1_lA(3,1) = (   -Az(i+2,j,k) + 8*Az(i+1,j,k)               &
                    - 8*Az(i-1,j,k) +   Az(i-2,j,k) ) / dx12

      d1_lA(1,2) = (   -Ax(i,j+2,k) + 8*Ax(i,j+1,k)               &
                    - 8*Ax(i,j-1,k) +   Ax(i,j-2,k) ) / dy12
      d1_lA(2,2) = (   -Ay(i,j+2,k) + 8*Ay(i,j+1,k)               &
                    - 8*Ay(i,j-1,k) +   Ay(i,j-2,k) ) / dy12
      d1_lA(3,2) = (   -Az(i,j+2,k) + 8*Az(i,j+1,k)               &
                    - 8*Az(i,j-1,k) +   Az(i,j-2,k) ) / dy12

      d1_lA(1,3) = (   -Ax(i,j,k+2) + 8*Ax(i,j,k+1)               &
                    - 8*Ax(i,j,k-1) +   Ax(i,j,k-2) ) / dz12
      d1_lA(2,3) = (   -Ay(i,j,k+2) + 8*Ay(i,j,k+1)               &
                    - 8*Ay(i,j,k-1) +   Ay(i,j,k-2) ) / dz12
      d1_lA(3,3) = (   -Az(i,j,k+2) + 8*Az(i,j,k+1)               &
                    - 8*Az(i,j,k-1) +   Az(i,j,k-2) ) / dz12

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

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

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

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

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

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

      


      

      
      d2_lA(1,1,1) = (  -Ax(i+2,j,k) + 16*Ax(i+1,j,k) - 30*Ax(i,j,k) &
                    + 16*Ax(i-1,j,k) -    Ax(i-2,j,k) ) / dxsq12
      d2_lA(2,1,1) = (  -Ay(i+2,j,k) + 16*Ay(i+1,j,k) - 30*Ay(i,j,k) &
                    + 16*Ay(i-1,j,k) -    Ay(i-2,j,k) ) / dxsq12
      d2_lA(3,1,1) = (  -Az(i+2,j,k) + 16*Az(i+1,j,k) - 30*Az(i,j,k) &
                    + 16*Az(i-1,j,k) -    Az(i-2,j,k) ) / dxsq12

      d2_lA(1,2,2) = (  -Ax(i,j+2,k) + 16*Ax(i,j+1,k) - 30*Ax(i,j,k) &
                    + 16*Ax(i,j-1,k) -    Ax(i,j-2,k) ) / dysq12
      d2_lA(2,2,2) = (  -Ay(i,j+2,k) + 16*Ay(i,j+1,k) - 30*Ay(i,j,k) &
                    + 16*Ay(i,j-1,k) -    Ay(i,j-2,k) ) / dysq12
      d2_lA(3,2,2) = (  -Az(i,j+2,k) + 16*Az(i,j+1,k) - 30*Az(i,j,k) &
                    + 16*Az(i,j-1,k) -    Az(i,j-2,k) ) / dysq12

      d2_lA(1,3,3) = (  -Ax(i,j,k+2) + 16*Ax(i,j,k+1) - 30*Ax(i,j,k) &
                    + 16*Ax(i,j,k-1) -    Ax(i,j,k-2) ) / dzsq12
      d2_lA(2,3,3) = (  -Ay(i,j,k+2) + 16*Ay(i,j,k+1) - 30*Ay(i,j,k) &
                    + 16*Ay(i,j,k-1) -    Ay(i,j,k-2) ) / dzsq12
      d2_lA(3,3,3) = (  -Az(i,j,k+2) + 16*Az(i,j,k+1) - 30*Az(i,j,k) &
                    + 16*Az(i,j,k-1) -    Az(i,j,k-2) ) / dzsq12

      d2_lA(1,1,2) = (  -Ax(i-2,j+2,k) +  8*Ax(i-1,j+2,k) -  8*Ax(i+1,j+2,k) +   Ax(i+2,j+2,k) &
                     + 8*Ax(i-2,j+1,k) - 64*Ax(i-1,j+1,k) + 64*Ax(i+1,j+1,k) - 8*Ax(i+2,j+1,k) &
                     - 8*Ax(i-2,j-1,k) + 64*Ax(i-1,j-1,k) - 64*Ax(i+1,j-1,k) + 8*Ax(i+2,j-1,k) &
                     +   Ax(i-2,j-2,k) -  8*Ax(i-1,j-2,k) +  8*Ax(i+1,j-2,k) -   Ax(i+2,j-2,k) ) / dxdy144
      d2_lA(2,1,2) = (  -Ay(i-2,j+2,k) +  8*Ay(i-1,j+2,k) -  8*Ay(i+1,j+2,k) +   Ay(i+2,j+2,k) &
                     + 8*Ay(i-2,j+1,k) - 64*Ay(i-1,j+1,k) + 64*Ay(i+1,j+1,k) - 8*Ay(i+2,j+1,k) &
                     - 8*Ay(i-2,j-1,k) + 64*Ay(i-1,j-1,k) - 64*Ay(i+1,j-1,k) + 8*Ay(i+2,j-1,k) &
                     +   Ay(i-2,j-2,k) -  8*Ay(i-1,j-2,k) +  8*Ay(i+1,j-2,k) -   Ay(i+2,j-2,k) ) / dxdy144
      d2_lA(3,1,2) = (  -Az(i-2,j+2,k) +  8*Az(i-1,j+2,k) -  8*Az(i+1,j+2,k) +   Az(i+2,j+2,k) &
                     + 8*Az(i-2,j+1,k) - 64*Az(i-1,j+1,k) + 64*Az(i+1,j+1,k) - 8*Az(i+2,j+1,k) &
                     - 8*Az(i-2,j-1,k) + 64*Az(i-1,j-1,k) - 64*Az(i+1,j-1,k) + 8*Az(i+2,j-1,k) &
                     +   Az(i-2,j-2,k) -  8*Az(i-1,j-2,k) +  8*Az(i+1,j-2,k) -   Az(i+2,j-2,k) ) / dxdy144

      d2_lA(1,1,3) = (  -Ax(i-2,j,k+2) +  8*Ax(i-1,j,k+2) -  8*Ax(i+1,j,k+2) +   Ax(i+2,j,k+2) &
                     + 8*Ax(i-2,j,k+1) - 64*Ax(i-1,j,k+1) + 64*Ax(i+1,j,k+1) - 8*Ax(i+2,j,k+1) &
                     - 8*Ax(i-2,j,k-1) + 64*Ax(i-1,j,k-1) - 64*Ax(i+1,j,k-1) + 8*Ax(i+2,j,k-1) &
                     +   Ax(i-2,j,k-2) -  8*Ax(i-1,j,k-2) +  8*Ax(i+1,j,k-2) -   Ax(i+2,j,k-2) ) / dxdz144
      d2_lA(2,1,3) = (  -Ay(i-2,j,k+2) +  8*Ay(i-1,j,k+2) -  8*Ay(i+1,j,k+2) +   Ay(i+2,j,k+2) &
                     + 8*Ay(i-2,j,k+1) - 64*Ay(i-1,j,k+1) + 64*Ay(i+1,j,k+1) - 8*Ay(i+2,j,k+1) &
                     - 8*Ay(i-2,j,k-1) + 64*Ay(i-1,j,k-1) - 64*Ay(i+1,j,k-1) + 8*Ay(i+2,j,k-1) &
                     +   Ay(i-2,j,k-2) -  8*Ay(i-1,j,k-2) +  8*Ay(i+1,j,k-2) -   Ay(i+2,j,k-2) ) / dxdz144
      d2_lA(3,1,3) = (  -Az(i-2,j,k+2) +  8*Az(i-1,j,k+2) -  8*Az(i+1,j,k+2) +   Az(i+2,j,k+2) &
                     + 8*Az(i-2,j,k+1) - 64*Az(i-1,j,k+1) + 64*Az(i+1,j,k+1) - 8*Az(i+2,j,k+1) &
                     - 8*Az(i-2,j,k-1) + 64*Az(i-1,j,k-1) - 64*Az(i+1,j,k-1) + 8*Az(i+2,j,k-1) &
                     +   Az(i-2,j,k-2) -  8*Az(i-1,j,k-2) +  8*Az(i+1,j,k-2) -   Az(i+2,j,k-2) ) / dxdz144

      d2_lA(1,2,3) = (  -Ax(i,j-2,k+2) +  8*Ax(i,j-1,k+2) -  8*Ax(i,j+1,k+2) +   Ax(i,j+2,k+2) &
                     + 8*Ax(i,j-2,k+1) - 64*Ax(i,j-1,k+1) + 64*Ax(i,j+1,k+1) - 8*Ax(i,j+2,k+1) &
                     - 8*Ax(i,j-2,k-1) + 64*Ax(i,j-1,k-1) - 64*Ax(i,j+1,k-1) + 8*Ax(i,j+2,k-1) &
                     +   Ax(i,j-2,k-2) -  8*Ax(i,j-1,k-2) +  8*Ax(i,j+1,k-2) -   Ax(i,j+2,k-2) ) / dydz144
      d2_lA(2,2,3) = (  -Ay(i,j-2,k+2) +  8*Ay(i,j-1,k+2) -  8*Ay(i,j+1,k+2) +   Ay(i,j+2,k+2) &
                     + 8*Ay(i,j-2,k+1) - 64*Ay(i,j-1,k+1) + 64*Ay(i,j+1,k+1) - 8*Ay(i,j+2,k+1) &
                     - 8*Ay(i,j-2,k-1) + 64*Ay(i,j-1,k-1) - 64*Ay(i,j+1,k-1) + 8*Ay(i,j+2,k-1) &
                     +   Ay(i,j-2,k-2) -  8*Ay(i,j-1,k-2) +  8*Ay(i,j+1,k-2) -   Ay(i,j+2,k-2) ) / dydz144
      d2_lA(3,2,3) = (  -Az(i,j-2,k+2) +  8*Az(i,j-1,k+2) -  8*Az(i,j+1,k+2) +   Az(i,j+2,k+2) &
                     + 8*Az(i,j-2,k+1) - 64*Az(i,j-1,k+1) + 64*Az(i,j+1,k+1) - 8*Az(i,j+2,k+1) &
                     - 8*Az(i,j-2,k-1) + 64*Az(i,j-1,k-1) - 64*Az(i,j+1,k-1) + 8*Az(i,j+2,k-1) &
                     +   Az(i,j-2,k-2) -  8*Az(i,j-1,k-2) +  8*Az(i,j+1,k-2) -   Az(i,j+2,k-2) ) / dydz144

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


      
      if( use_advection_stencils /= 0 ) then

        di = int( sign( one, beta(1) ) )
        dj = int( sign( one, beta(2) ) )
        dk = int( sign( one, beta(3) ) )

        
        d1_f(1) = di * ( -3*Ex(i-di,j,k) - 10*Ex(i,j,k) + 18*Ex(i+di,j,k) &
                        - 6*Ex(i+2*di,j,k) + Ex(i+3*di,j,k)) / dx12
        d1_f(2) = dj * ( -3*Ex(i,j-dj,k) - 10*Ex(i,j,k) + 18*Ex(i,j+dj,k) &
                        - 6*Ex(i,j+2*dj,k) + Ex(i,j+3*dj,k)) / dy12
        d1_f(3) = dk * ( -3*Ex(i,j,k-dk) - 10*Ex(i,j,k) + 18*Ex(i,j,k+dk) &
                        - 6*Ex(i,j,k+2*dk) + Ex(i,j,k+3*dk)) / dz12
        ad1_lE(1) = beta(1)*d1_f(1) + beta(2)*d1_f(2) + beta(3)*d1_f(3)

        d1_f(1) = di * ( -3*Ey(i-di,j,k) - 10*Ey(i,j,k) + 18*Ey(i+di,j,k) &
                        - 6*Ey(i+2*di,j,k) + Ey(i+3*di,j,k)) / dx12
        d1_f(2) = dj * ( -3*Ey(i,j-dj,k) - 10*Ey(i,j,k) + 18*Ey(i,j+dj,k) &
                        - 6*Ey(i,j+2*dj,k) + Ey(i,j+3*dj,k)) / dy12
        d1_f(3) = dk * ( -3*Ey(i,j,k-dk) - 10*Ey(i,j,k) + 18*Ey(i,j,k+dk) &
                        - 6*Ey(i,j,k+2*dk) + Ey(i,j,k+3*dk)) / dz12
        ad1_lE(2) = beta(1)*d1_f(1) + beta(2)*d1_f(2) + beta(3)*d1_f(3)

        d1_f(1) = di * ( -3*Ez(i-di,j,k) - 10*Ez(i,j,k) + 18*Ez(i+di,j,k) &
                        - 6*Ez(i+2*di,j,k) + Ez(i+3*di,j,k)) / dx12
        d1_f(2) = dj * ( -3*Ez(i,j-dj,k) - 10*Ez(i,j,k) + 18*Ez(i,j+dj,k) &
                        - 6*Ez(i,j+2*dj,k) + Ez(i,j+3*dj,k)) / dy12
        d1_f(3) = dk * ( -3*Ez(i,j,k-dk) - 10*Ez(i,j,k) + 18*Ez(i,j,k+dk) &
                        - 6*Ez(i,j,k+2*dk) + Ez(i,j,k+3*dk)) / dz12
        ad1_lE(3) = beta(1)*d1_f(1) + beta(2)*d1_f(2) + beta(3)*d1_f(3)

        
        d1_f(1) = di * ( -3*Ax(i-di,j,k) - 10*Ax(i,j,k) + 18*Ax(i+di,j,k) &
                        - 6*Ax(i+2*di,j,k) + Ax(i+3*di,j,k)) / dx12
        d1_f(2) = dj * ( -3*Ax(i,j-dj,k) - 10*Ax(i,j,k) + 18*Ax(i,j+dj,k) &
                        - 6*Ax(i,j+2*dj,k) + Ax(i,j+3*dj,k)) / dy12
        d1_f(3) = dk * ( -3*Ax(i,j,k-dk) - 10*Ax(i,j,k) + 18*Ax(i,j,k+dk) &
                        - 6*Ax(i,j,k+2*dk) + Ax(i,j,k+3*dk)) / dz12
        ad1_lA(1) = beta(1)*d1_f(1) + beta(2)*d1_f(2) + beta(3)*d1_f(3)

        d1_f(1) = di * ( -3*Ay(i-di,j,k) - 10*Ay(i,j,k) + 18*Ay(i+di,j,k) &
                        - 6*Ay(i+2*di,j,k) + Ay(i+3*di,j,k)) / dx12
        d1_f(2) = dj * ( -3*Ay(i,j-dj,k) - 10*Ay(i,j,k) + 18*Ay(i,j+dj,k) &
                        - 6*Ay(i,j+2*dj,k) + Ay(i,j+3*dj,k)) / dy12
        d1_f(3) = dk * ( -3*Ay(i,j,k-dk) - 10*Ay(i,j,k) + 18*Ay(i,j,k+dk) &
                        - 6*Ay(i,j,k+2*dk) + Ay(i,j,k+3*dk)) / dz12
        ad1_lA(2) = beta(1)*d1_f(1) + beta(2)*d1_f(2) + beta(3)*d1_f(3)

        d1_f(1) = di * ( -3*Az(i-di,j,k) - 10*Az(i,j,k) + 18*Az(i+di,j,k) &
                        - 6*Az(i+2*di,j,k) + Az(i+3*di,j,k)) / dx12
        d1_f(2) = dj * ( -3*Az(i,j-dj,k) - 10*Az(i,j,k) + 18*Az(i,j+dj,k) &
                        - 6*Az(i,j+2*dj,k) + Az(i,j+3*dj,k)) / dy12
        d1_f(3) = dk * ( -3*Az(i,j,k-dk) - 10*Az(i,j,k) + 18*Az(i,j,k+dk) &
                        - 6*Az(i,j,k+2*dk) + Az(i,j,k+3*dk)) / dz12
        ad1_lA(3) = beta(1)*d1_f(1) + beta(2)*d1_f(2) + beta(3)*d1_f(3)

        
        d1_f(1) = di * ( -3*Zeta(i-di,j,k) - 10*Zeta(i,j,k) + 18*Zeta(i+di,j,k)  &
                        - 6*Zeta(i+2*di,j,k) + Zeta(i+3*di,j,k)) / dx12
        d1_f(2) = dj * ( -3*Zeta(i,j-dj,k) - 10*Zeta(i,j,k) + 18*Zeta(i,j+dj,k)  &
                        - 6*Zeta(i,j+2*dj,k) + Zeta(i,j+3*dj,k)) / dy12
        d1_f(3) = dk * ( -3*Zeta(i,j,k-dk) - 10*Zeta(i,j,k) + 18*Zeta(i,j,k+dk)  &
                        - 6*Zeta(i,j,k+2*dk) + Zeta(i,j,k+3*dk)) / dz12
        ad1_lZeta = beta(1)*d1_f(1) + beta(2)*d1_f(2) + beta(3)*d1_f(3)

        
        d1_f(1) = di * ( -3*Aphi(i-di,j,k) - 10*Aphi(i,j,k) + 18*Aphi(i+di,j,k)  &
                        - 6*Aphi(i+2*di,j,k) + Aphi(i+3*di,j,k)) / dx12
        d1_f(2) = dj * ( -3*Aphi(i,j-dj,k) - 10*Aphi(i,j,k) + 18*Aphi(i,j+dj,k)  &
                        - 6*Aphi(i,j+2*dj,k) + Aphi(i,j+3*dj,k)) / dy12
        d1_f(3) = dk * ( -3*Aphi(i,j,k-dk) - 10*Aphi(i,j,k) + 18*Aphi(i,j,k+dk)  &
                        - 6*Aphi(i,j,k+2*dk) + Aphi(i,j,k+3*dk)) / dz12
        ad1_lAphi = beta(1)*d1_f(1) + beta(2)*d1_f(2) + beta(3)*d1_f(3)

      else

        
        ad1_lE(1) = beta(1)*d1_lE(1,1) + beta(2)*d1_lE(1,2) + beta(3)*d1_lE(1,3)
        ad1_lE(2) = beta(1)*d1_lE(2,1) + beta(2)*d1_lE(2,2) + beta(3)*d1_lE(2,3)
        ad1_lE(3) = beta(1)*d1_lE(3,1) + beta(2)*d1_lE(3,2) + beta(3)*d1_lE(3,3)

        
        ad1_lA(1) = beta(1)*d1_lA(1,1) + beta(2)*d1_lA(1,2) + beta(3)*d1_lA(1,3)
        ad1_lA(2) = beta(1)*d1_lA(2,1) + beta(2)*d1_lA(2,2) + beta(3)*d1_lA(2,3)
        ad1_lA(3) = beta(1)*d1_lA(3,1) + beta(2)*d1_lA(3,2) + beta(3)*d1_lA(3,3)

        
        ad1_lZeta = beta(1)*d1_lZeta(1) + beta(2)*d1_lZeta(2) + beta(3)*d1_lZeta(3)

        
        ad1_lAphi = beta(1)*d1_lAphi(1) + beta(2)*d1_lAphi(2) + beta(3)*d1_lAphi(3)

      end if
      

    else if (derivs_order == 6) then

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

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

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

      
      d1_hh(1,1,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_hh(1,2,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_hh(1,3,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_hh(2,2,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_hh(2,3,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_hh(3,3,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_hh(1,1,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_hh(1,2,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_hh(1,3,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_hh(2,2,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_hh(2,3,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_hh(3,3,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_hh(1,1,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_hh(1,2,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_hh(1,3,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_hh(2,2,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_hh(2,3,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_hh(3,3,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_hh(2,1,:) = d1_hh(1,2,:)
      d1_hh(3,1,:) = d1_hh(1,3,:)
      d1_hh(3,2,:) = d1_hh(2,3,:)

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

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

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

      
      d1_beta(1,1) = (  betax(i+3,j,k) - 9*betax(i+2,j,k) + 45*betax(i+1,j,k) &
                      - betax(i-3,j,k) + 9*betax(i-2,j,k) - 45*betax(i-1,j,k) ) * odx60
      d1_beta(2,1) = (  betay(i+3,j,k) - 9*betay(i+2,j,k) + 45*betay(i+1,j,k) &
                      - betay(i-3,j,k) + 9*betay(i-2,j,k) - 45*betay(i-1,j,k) ) * odx60
      d1_beta(3,1) = (  betaz(i+3,j,k) - 9*betaz(i+2,j,k) + 45*betaz(i+1,j,k) &
                      - betaz(i-3,j,k) + 9*betaz(i-2,j,k) - 45*betaz(i-1,j,k) ) * odx60

      d1_beta(1,2) = (  betax(i,j+3,k) - 9*betax(i,j+2,k) + 45*betax(i,j+1,k) &
                      - betax(i,j-3,k) + 9*betax(i,j-2,k) - 45*betax(i,j-1,k) ) * ody60
      d1_beta(2,2) = (  betay(i,j+3,k) - 9*betay(i,j+2,k) + 45*betay(i,j+1,k) &
                      - betay(i,j-3,k) + 9*betay(i,j-2,k) - 45*betay(i,j-1,k) ) * ody60
      d1_beta(3,2) = (  betaz(i,j+3,k) - 9*betaz(i,j+2,k) + 45*betaz(i,j+1,k) &
                      - betaz(i,j-3,k) + 9*betaz(i,j-2,k) - 45*betaz(i,j-1,k) ) * ody60

      d1_beta(1,3) = (  betax(i,j,k+3) - 9*betax(i,j,k+2) + 45*betax(i,j,k+1) &
                      - betax(i,j,k-3) + 9*betax(i,j,k-2) - 45*betax(i,j,k-1) ) * odz60
      d1_beta(2,3) = (  betay(i,j,k+3) - 9*betay(i,j,k+2) + 45*betay(i,j,k+1) &
                      - betay(i,j,k-3) + 9*betay(i,j,k-2) - 45*betay(i,j,k-1) ) * odz60
      d1_beta(3,3) = (  betaz(i,j,k+3) - 9*betaz(i,j,k+2) + 45*betaz(i,j,k+1) &
                      - betaz(i,j,k-3) + 9*betaz(i,j,k-2) - 45*betaz(i,j,k-1) ) * odz60

      
      d1_lE(1,1) = (  Ex(i+3,j,k) - 9*Ex(i+2,j,k) + 45*Ex(i+1,j,k) &
                    - Ex(i-3,j,k) + 9*Ex(i-2,j,k) - 45*Ex(i-1,j,k) ) * odx60
      d1_lE(2,1) = (  Ey(i+3,j,k) - 9*Ey(i+2,j,k) + 45*Ey(i+1,j,k) &
                    - Ey(i-3,j,k) + 9*Ey(i-2,j,k) - 45*Ey(i-1,j,k) ) * odx60
      d1_lE(3,1) = (  Ez(i+3,j,k) - 9*Ez(i+2,j,k) + 45*Ez(i+1,j,k) &
                    - Ez(i-3,j,k) + 9*Ez(i-2,j,k) - 45*Ez(i-1,j,k) ) * odx60

      d1_lE(1,2) = (  Ex(i,j+3,k) - 9*Ex(i,j+2,k) + 45*Ex(i,j+1,k) &
                    - Ex(i,j-3,k) + 9*Ex(i,j-2,k) - 45*Ex(i,j-1,k) ) * ody60
      d1_lE(2,2) = (  Ey(i,j+3,k) - 9*Ey(i,j+2,k) + 45*Ey(i,j+1,k) &
                    - Ey(i,j-3,k) + 9*Ey(i,j-2,k) - 45*Ey(i,j-1,k) ) * ody60
      d1_lE(3,2) = (  Ez(i,j+3,k) - 9*Ez(i,j+2,k) + 45*Ez(i,j+1,k) &
                    - Ez(i,j-3,k) + 9*Ez(i,j-2,k) - 45*Ez(i,j-1,k) ) * ody60

      d1_lE(1,3) = (  Ex(i,j,k+3) - 9*Ex(i,j,k+2) + 45*Ex(i,j,k+1) &
                    - Ex(i,j,k-3) + 9*Ex(i,j,k-2) - 45*Ex(i,j,k-1) ) * odz60
      d1_lE(2,3) = (  Ey(i,j,k+3) - 9*Ey(i,j,k+2) + 45*Ey(i,j,k+1) &
                    - Ey(i,j,k-3) + 9*Ey(i,j,k-2) - 45*Ey(i,j,k-1) ) * odz60
      d1_lE(3,3) = (  Ez(i,j,k+3) - 9*Ez(i,j,k+2) + 45*Ez(i,j,k+1) &
                    - Ez(i,j,k-3) + 9*Ez(i,j,k-2) - 45*Ez(i,j,k-1) ) * odz60

      
      d1_lA(1,1) = (  Ax(i+3,j,k) - 9*Ax(i+2,j,k) + 45*Ax(i+1,j,k) &
                    - Ax(i-3,j,k) + 9*Ax(i-2,j,k) - 45*Ax(i-1,j,k) ) * odx60
      d1_lA(2,1) = (  Ay(i+3,j,k) - 9*Ay(i+2,j,k) + 45*Ay(i+1,j,k) &
                    - Ay(i-3,j,k) + 9*Ay(i-2,j,k) - 45*Ay(i-1,j,k) ) * odx60
      d1_lA(3,1) = (  Az(i+3,j,k) - 9*Az(i+2,j,k) + 45*Az(i+1,j,k) &
                    - Az(i-3,j,k) + 9*Az(i-2,j,k) - 45*Az(i-1,j,k) ) * odx60

      d1_lA(1,2) = (  Ax(i,j+3,k) - 9*Ax(i,j+2,k) + 45*Ax(i,j+1,k) &
                    - Ax(i,j-3,k) + 9*Ax(i,j-2,k) - 45*Ax(i,j-1,k) ) * ody60
      d1_lA(2,2) = (  Ay(i,j+3,k) - 9*Ay(i,j+2,k) + 45*Ay(i,j+1,k) &
                    - Ay(i,j-3,k) + 9*Ay(i,j-2,k) - 45*Ay(i,j-1,k) ) * ody60
      d1_lA(3,2) = (  Az(i,j+3,k) - 9*Az(i,j+2,k) + 45*Az(i,j+1,k) &
                    - Az(i,j-3,k) + 9*Az(i,j-2,k) - 45*Az(i,j-1,k) ) * ody60

      d1_lA(1,3) = (  Ax(i,j,k+3) - 9*Ax(i,j,k+2) + 45*Ax(i,j,k+1) &
                    - Ax(i,j,k-3) + 9*Ax(i,j,k-2) - 45*Ax(i,j,k-1) ) * odz60
      d1_lA(2,3) = (  Ay(i,j,k+3) - 9*Ay(i,j,k+2) + 45*Ay(i,j,k+1) &
                    - Ay(i,j,k-3) + 9*Ay(i,j,k-2) - 45*Ay(i,j,k-1) ) * odz60
      d1_lA(3,3) = (  Az(i,j,k+3) - 9*Az(i,j,k+2) + 45*Az(i,j,k+1) &
                    - Az(i,j,k-3) + 9*Az(i,j,k-2) - 45*Az(i,j,k-1) ) * odz60

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

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

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

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

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

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

      


      

      
      d2_lA(1,1,1) = (  2*Ax(i+3,j,k) - 27*Ax(i+2,j,k) + 270*Ax(i+1,j,k) - 490*Ax(i,j,k)&
                      + 2*Ax(i-3,j,k) - 27*Ax(i-2,j,k) + 270*Ax(i-1,j,k) ) * odxsq180
      d2_lA(2,1,1) = (  2*Ay(i+3,j,k) - 27*Ay(i+2,j,k) + 270*Ay(i+1,j,k) - 490*Ay(i,j,k)&
                      + 2*Ay(i-3,j,k) - 27*Ay(i-2,j,k) + 270*Ay(i-1,j,k) ) * odxsq180
      d2_lA(3,1,1) = (  2*Az(i+3,j,k) - 27*Az(i+2,j,k) + 270*Az(i+1,j,k) - 490*Az(i,j,k)&
                      + 2*Az(i-3,j,k) - 27*Az(i-2,j,k) + 270*Az(i-1,j,k) ) * odxsq180

      d2_lA(1,2,2) = (  2*Ax(i,j+3,k) - 27*Ax(i,j+2,k) + 270*Ax(i,j+1,k) - 490*Ax(i,j,k)&
                      + 2*Ax(i,j-3,k) - 27*Ax(i,j-2,k) + 270*Ax(i,j-1,k) ) * odysq180
      d2_lA(2,2,2) = (  2*Ay(i,j+3,k) - 27*Ay(i,j+2,k) + 270*Ay(i,j+1,k) - 490*Ay(i,j,k)&
                      + 2*Ay(i,j-3,k) - 27*Ay(i,j-2,k) + 270*Ay(i,j-1,k) ) * odysq180
      d2_lA(3,2,2) = (  2*Az(i,j+3,k) - 27*Az(i,j+2,k) + 270*Az(i,j+1,k) - 490*Az(i,j,k)&
                      + 2*Az(i,j-3,k) - 27*Az(i,j-2,k) + 270*Az(i,j-1,k) ) * odysq180

      d2_lA(1,3,3) = (  2*Ax(i,j,k+3) - 27*Ax(i,j,k+2) + 270*Ax(i,j,k+1) - 490*Ax(i,j,k)&
                      + 2*Ax(i,j,k-3) - 27*Ax(i,j,k-2) + 270*Ax(i,j,k-1) ) * odzsq180
      d2_lA(2,3,3) = (  2*Ay(i,j,k+3) - 27*Ay(i,j,k+2) + 270*Ay(i,j,k+1) - 490*Ay(i,j,k)&
                      + 2*Ay(i,j,k-3) - 27*Ay(i,j,k-2) + 270*Ay(i,j,k-1) ) * odzsq180
      d2_lA(3,3,3) = (  2*Az(i,j,k+3) - 27*Az(i,j,k+2) + 270*Az(i,j,k+1) - 490*Az(i,j,k)&
                      + 2*Az(i,j,k-3) - 27*Az(i,j,k-2) + 270*Az(i,j,k-1) ) * odzsq180

      d2_lA(1,1,2) = (    -Ax(i-3,j+3,k) +   9*Ax(i-2,j+3,k) -   45*Ax(i-1,j+3,k) +   45*Ax(i+1,j+3,k) -   9*Ax(i+2,j+3,k) +    Ax(&
  &i+3,j+3,k) &
                      +  9*Ax(i-3,j+2,k) -  81*Ax(i-2,j+2,k) +  405*Ax(i-1,j+2,k) -  405*Ax(i+1,j+2,k) +  81*Ax(i+2,j+2,k) -  9*Ax(&
  &i+3,j+2,k) &
                      - 45*Ax(i-3,j+1,k) + 405*Ax(i-2,j+1,k) - 2025*Ax(i-1,j+1,k) + 2025*Ax(i+1,j+1,k) - 405*Ax(i+2,j+1,k) + 45*Ax(&
  &i+3,j+1,k) &
                      + 45*Ax(i-3,j-1,k) - 405*Ax(i-2,j-1,k) + 2025*Ax(i-1,j-1,k) - 2025*Ax(i+1,j-1,k) + 405*Ax(i+2,j-1,k) - 45*Ax(&
  &i+3,j-1,k) &
                      -  9*Ax(i-3,j-2,k) +  81*Ax(i-2,j-2,k) -  405*Ax(i-1,j-2,k) +  405*Ax(i+1,j-2,k) -  81*Ax(i+2,j-2,k) +  9*Ax(&
  &i+3,j-2,k) &
                      +    Ax(i-3,j-3,k) -   9*Ax(i-2,j-3,k) +   45*Ax(i-1,j-3,k) -   45*Ax(i+1,j-3,k) +   9*Ax(i+2,j-3,k) -    Ax(&
  &i+3,j-3,k) ) * odxdy3600
      d2_lA(2,1,2) = (    -Ay(i-3,j+3,k) +   9*Ay(i-2,j+3,k) -   45*Ay(i-1,j+3,k) +   45*Ay(i+1,j+3,k) -   9*Ay(i+2,j+3,k) +    Ay(&
  &i+3,j+3,k) &
                      +  9*Ay(i-3,j+2,k) -  81*Ay(i-2,j+2,k) +  405*Ay(i-1,j+2,k) -  405*Ay(i+1,j+2,k) +  81*Ay(i+2,j+2,k) -  9*Ay(&
  &i+3,j+2,k) &
                      - 45*Ay(i-3,j+1,k) + 405*Ay(i-2,j+1,k) - 2025*Ay(i-1,j+1,k) + 2025*Ay(i+1,j+1,k) - 405*Ay(i+2,j+1,k) + 45*Ay(&
  &i+3,j+1,k) &
                      + 45*Ay(i-3,j-1,k) - 405*Ay(i-2,j-1,k) + 2025*Ay(i-1,j-1,k) - 2025*Ay(i+1,j-1,k) + 405*Ay(i+2,j-1,k) - 45*Ay(&
  &i+3,j-1,k) &
                      -  9*Ay(i-3,j-2,k) +  81*Ay(i-2,j-2,k) -  405*Ay(i-1,j-2,k) +  405*Ay(i+1,j-2,k) -  81*Ay(i+2,j-2,k) +  9*Ay(&
  &i+3,j-2,k) &
                      +    Ay(i-3,j-3,k) -   9*Ay(i-2,j-3,k) +   45*Ay(i-1,j-3,k) -   45*Ay(i+1,j-3,k) +   9*Ay(i+2,j-3,k) -    Ay(&
  &i+3,j-3,k) ) * odxdy3600
      d2_lA(3,1,2) = (    -Az(i-3,j+3,k) +   9*Az(i-2,j+3,k) -   45*Az(i-1,j+3,k) +   45*Az(i+1,j+3,k) -   9*Az(i+2,j+3,k) +    Az(&
  &i+3,j+3,k) &
                      +  9*Az(i-3,j+2,k) -  81*Az(i-2,j+2,k) +  405*Az(i-1,j+2,k) -  405*Az(i+1,j+2,k) +  81*Az(i+2,j+2,k) -  9*Az(&
  &i+3,j+2,k) &
                      - 45*Az(i-3,j+1,k) + 405*Az(i-2,j+1,k) - 2025*Az(i-1,j+1,k) + 2025*Az(i+1,j+1,k) - 405*Az(i+2,j+1,k) + 45*Az(&
  &i+3,j+1,k) &
                      + 45*Az(i-3,j-1,k) - 405*Az(i-2,j-1,k) + 2025*Az(i-1,j-1,k) - 2025*Az(i+1,j-1,k) + 405*Az(i+2,j-1,k) - 45*Az(&
  &i+3,j-1,k) &
                      -  9*Az(i-3,j-2,k) +  81*Az(i-2,j-2,k) -  405*Az(i-1,j-2,k) +  405*Az(i+1,j-2,k) -  81*Az(i+2,j-2,k) +  9*Az(&
  &i+3,j-2,k) &
                      +    Az(i-3,j-3,k) -   9*Az(i-2,j-3,k) +   45*Az(i-1,j-3,k) -   45*Az(i+1,j-3,k) +   9*Az(i+2,j-3,k) -    Az(&
  &i+3,j-3,k) ) * odxdy3600


      d2_lA(1,1,3) = (    -Ax(i-3,j,k+3) +   9*Ax(i-2,j,k+3) -   45*Ax(i-1,j,k+3) +   45*Ax(i+1,j,k+3) -   9*Ax(i+2,j,k+3) +    Ax(&
  &i+3,j,k+3) &
                      +  9*Ax(i-3,j,k+2) -  81*Ax(i-2,j,k+2) +  405*Ax(i-1,j,k+2) -  405*Ax(i+1,j,k+2) +  81*Ax(i+2,j,k+2) -  9*Ax(&
  &i+3,j,k+2) &
                      - 45*Ax(i-3,j,k+1) + 405*Ax(i-2,j,k+1) - 2025*Ax(i-1,j,k+1) + 2025*Ax(i+1,j,k+1) - 405*Ax(i+2,j,k+1) + 45*Ax(&
  &i+3,j,k+1) &
                      + 45*Ax(i-3,j,k-1) - 405*Ax(i-2,j,k-1) + 2025*Ax(i-1,j,k-1) - 2025*Ax(i+1,j,k-1) + 405*Ax(i+2,j,k-1) - 45*Ax(&
  &i+3,j,k-1) &
                      -  9*Ax(i-3,j,k-2) +  81*Ax(i-2,j,k-2) -  405*Ax(i-1,j,k-2) +  405*Ax(i+1,j,k-2) -  81*Ax(i+2,j,k-2) +  9*Ax(&
  &i+3,j,k-2) &
                      +    Ax(i-3,j,k-3) -   9*Ax(i-2,j,k-3) +   45*Ax(i-1,j,k-3) -   45*Ax(i+1,j,k-3) +   9*Ax(i+2,j,k-3) -    Ax(&
  &i+3,j,k-3) ) * odxdz3600
      d2_lA(2,1,3) = (    -Ay(i-3,j,k+3) +   9*Ay(i-2,j,k+3) -   45*Ay(i-1,j,k+3) +   45*Ay(i+1,j,k+3) -   9*Ay(i+2,j,k+3) +    Ay(&
  &i+3,j,k+3) &
                      +  9*Ay(i-3,j,k+2) -  81*Ay(i-2,j,k+2) +  405*Ay(i-1,j,k+2) -  405*Ay(i+1,j,k+2) +  81*Ay(i+2,j,k+2) -  9*Ay(&
  &i+3,j,k+2) &
                      - 45*Ay(i-3,j,k+1) + 405*Ay(i-2,j,k+1) - 2025*Ay(i-1,j,k+1) + 2025*Ay(i+1,j,k+1) - 405*Ay(i+2,j,k+1) + 45*Ay(&
  &i+3,j,k+1) &
                      + 45*Ay(i-3,j,k-1) - 405*Ay(i-2,j,k-1) + 2025*Ay(i-1,j,k-1) - 2025*Ay(i+1,j,k-1) + 405*Ay(i+2,j,k-1) - 45*Ay(&
  &i+3,j,k-1) &
                      -  9*Ay(i-3,j,k-2) +  81*Ay(i-2,j,k-2) -  405*Ay(i-1,j,k-2) +  405*Ay(i+1,j,k-2) -  81*Ay(i+2,j,k-2) +  9*Ay(&
  &i+3,j,k-2) &
                      +    Ay(i-3,j,k-3) -   9*Ay(i-2,j,k-3) +   45*Ay(i-1,j,k-3) -   45*Ay(i+1,j,k-3) +   9*Ay(i+2,j,k-3) -    Ay(&
  &i+3,j,k-3) ) * odxdz3600
      d2_lA(3,1,3) = (    -Az(i-3,j,k+3) +   9*Az(i-2,j,k+3) -   45*Az(i-1,j,k+3) +   45*Az(i+1,j,k+3) -   9*Az(i+2,j,k+3) +    Az(&
  &i+3,j,k+3) &
                      +  9*Az(i-3,j,k+2) -  81*Az(i-2,j,k+2) +  405*Az(i-1,j,k+2) -  405*Az(i+1,j,k+2) +  81*Az(i+2,j,k+2) -  9*Az(&
  &i+3,j,k+2) &
                      - 45*Az(i-3,j,k+1) + 405*Az(i-2,j,k+1) - 2025*Az(i-1,j,k+1) + 2025*Az(i+1,j,k+1) - 405*Az(i+2,j,k+1) + 45*Az(&
  &i+3,j,k+1) &
                      + 45*Az(i-3,j,k-1) - 405*Az(i-2,j,k-1) + 2025*Az(i-1,j,k-1) - 2025*Az(i+1,j,k-1) + 405*Az(i+2,j,k-1) - 45*Az(&
  &i+3,j,k-1) &
                      -  9*Az(i-3,j,k-2) +  81*Az(i-2,j,k-2) -  405*Az(i-1,j,k-2) +  405*Az(i+1,j,k-2) -  81*Az(i+2,j,k-2) +  9*Az(&
  &i+3,j,k-2) &
                      +    Az(i-3,j,k-3) -   9*Az(i-2,j,k-3) +   45*Az(i-1,j,k-3) -   45*Az(i+1,j,k-3) +   9*Az(i+2,j,k-3) -    Az(&
  &i+3,j,k-3) ) * odxdz3600

      d2_lA(1,2,3) = (    -Ax(i,j-3,k+3) +   9*Ax(i,j-2,k+3) -   45*Ax(i,j-1,k+3) +   45*Ax(i,j+1,k+3) -   9*Ax(i,j+2,k+3) +    Ax(&
  &i,j+3,k+3) &
                      +  9*Ax(i,j-3,k+2) -  81*Ax(i,j-2,k+2) +  405*Ax(i,j-1,k+2) -  405*Ax(i,j+1,k+2) +  81*Ax(i,j+2,k+2) -  9*Ax(&
  &i,j+3,k+2) &
                      - 45*Ax(i,j-3,k+1) + 405*Ax(i,j-2,k+1) - 2025*Ax(i,j-1,k+1) + 2025*Ax(i,j+1,k+1) - 405*Ax(i,j+2,k+1) + 45*Ax(&
  &i,j+3,k+1) &
                      + 45*Ax(i,j-3,k-1) - 405*Ax(i,j-2,k-1) + 2025*Ax(i,j-1,k-1) - 2025*Ax(i,j+1,k-1) + 405*Ax(i,j+2,k-1) - 45*Ax(&
  &i,j+3,k-1) &
                      -  9*Ax(i,j-3,k-2) +  81*Ax(i,j-2,k-2) -  405*Ax(i,j-1,k-2) +  405*Ax(i,j+1,k-2) -  81*Ax(i,j+2,k-2) +  9*Ax(&
  &i,j+3,k-2) &
                      +    Ax(i,j-3,k-3) -   9*Ax(i,j-2,k-3) +   45*Ax(i,j-1,k-3) -   45*Ax(i,j+1,k-3) +   9*Ax(i,j+2,k-3) -    Ax(&
  &i,j+3,k-3) ) * odydz3600
      d2_lA(2,2,3) = (    -Ay(i,j-3,k+3) +   9*Ay(i,j-2,k+3) -   45*Ay(i,j-1,k+3) +   45*Ay(i,j+1,k+3) -   9*Ay(i,j+2,k+3) +    Ay(&
  &i,j+3,k+3) &
                      +  9*Ay(i,j-3,k+2) -  81*Ay(i,j-2,k+2) +  405*Ay(i,j-1,k+2) -  405*Ay(i,j+1,k+2) +  81*Ay(i,j+2,k+2) -  9*Ay(&
  &i,j+3,k+2) &
                      - 45*Ay(i,j-3,k+1) + 405*Ay(i,j-2,k+1) - 2025*Ay(i,j-1,k+1) + 2025*Ay(i,j+1,k+1) - 405*Ay(i,j+2,k+1) + 45*Ay(&
  &i,j+3,k+1) &
                      + 45*Ay(i,j-3,k-1) - 405*Ay(i,j-2,k-1) + 2025*Ay(i,j-1,k-1) - 2025*Ay(i,j+1,k-1) + 405*Ay(i,j+2,k-1) - 45*Ay(&
  &i,j+3,k-1) &
                      -  9*Ay(i,j-3,k-2) +  81*Ay(i,j-2,k-2) -  405*Ay(i,j-1,k-2) +  405*Ay(i,j+1,k-2) -  81*Ay(i,j+2,k-2) +  9*Ay(&
  &i,j+3,k-2) &
                     +    Ay(i,j-3,k-3) -   9*Ay(i,j-2,k-3) +   45*Ay(i,j-1,k-3) -   45*Ay(i,j+1,k-3) +   9*Ay(i,j+2,k-3) -    Ay(i&
  &,j+3,k-3) ) * odydz3600
      d2_lA(3,2,3) = (    -Az(i,j-3,k+3) +   9*Az(i,j-2,k+3) -   45*Az(i,j-1,k+3) +   45*Az(i,j+1,k+3) -   9*Az(i,j+2,k+3) +    Az(&
  &i,j+3,k+3) &
                      +  9*Az(i,j-3,k+2) -  81*Az(i,j-2,k+2) +  405*Az(i,j-1,k+2) -  405*Az(i,j+1,k+2) +  81*Az(i,j+2,k+2) -  9*Az(&
  &i,j+3,k+2) &
                      - 45*Az(i,j-3,k+1) + 405*Az(i,j-2,k+1) - 2025*Az(i,j-1,k+1) + 2025*Az(i,j+1,k+1) - 405*Az(i,j+2,k+1) + 45*Az(&
  &i,j+3,k+1) &
                      + 45*Az(i,j-3,k-1) - 405*Az(i,j-2,k-1) + 2025*Az(i,j-1,k-1) - 2025*Az(i,j+1,k-1) + 405*Az(i,j+2,k-1) - 45*Az(&
  &i,j+3,k-1) &
                      -  9*Az(i,j-3,k-2) +  81*Az(i,j-2,k-2) -  405*Az(i,j-1,k-2) +  405*Az(i,j+1,k-2) -  81*Az(i,j+2,k-2) +  9*Az(&
  &i,j+3,k-2) &
                      +    Az(i,j-3,k-3) -   9*Az(i,j-2,k-3) +   45*Az(i,j-1,k-3) -   45*Az(i,j+1,k-3) +   9*Az(i,j+2,k-3) -    Az(&
  &i,j+3,k-3) ) * odydz3600

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


      
      if( use_advection_stencils /= 0 ) then

        di = int( sign( one, beta(1) ) )
        dj = int( sign( one, beta(2) ) )
        dk = int( sign( one, beta(3) ) )

        
        d1_f(1) = di * (   2*Ex(i-2*di,j,k) - 24*Ex(i-di,j,k) - 35*Ex(i,j,k) + 80*Ex(i+di,j,k) &
                        - 30*Ex(i+2*di,j,k) + 8*Ex(i+3*di,j,k) - Ex(i+4*di,j,k) ) * odx60
        d1_f(2) = dj * (   2*Ex(i,j-2*dj,k) - 24*Ex(i,j-dj,k) - 35*Ex(i,j,k) + 80*Ex(i,j+dj,k) &
                        - 30*Ex(i,j+2*dj,k) + 8*Ex(i,j+3*dj,k) - Ex(i,j+4*dj,k) ) * ody60
        d1_f(3) = dk * (   2*Ex(i,j,k-2*dk) - 24*Ex(i,j,k-dk) - 35*Ex(i,j,k) + 80*Ex(i,j,k+dk) &
                        - 30*Ex(i,j,k+2*dk) + 8*Ex(i,j,k+3*dk) - Ex(i,j,k+4*dk) ) * odz60
        ad1_lE(1) = beta(1)*d1_f(1) + beta(2)*d1_f(2) + beta(3)*d1_f(3)

        d1_f(1) = di * (   2*Ey(i-2*di,j,k) - 24*Ey(i-di,j,k) - 35*Ey(i,j,k) + 80*Ey(i+di,j,k) &
                        - 30*Ey(i+2*di,j,k) +  8*Ey(i+3*di,j,k) -  Ey(i+4*di,j,k) ) * odx60
        d1_f(2) = dj * (   2*Ey(i,j-2*dj,k) - 24*Ey(i,j-dj,k) - 35*Ey(i,j,k) + 80*Ey(i,j+dj,k) &
                        - 30*Ey(i,j+2*dj,k) +  8*Ey(i,j+3*dj,k) -  Ey(i,j+4*dj,k) ) * ody60
        d1_f(3) = dk * (   2*Ey(i,j,k-2*dk) - 24*Ey(i,j,k-dk) - 35*Ey(i,j,k) + 80*Ey(i,j,k+dk) &
                        - 30*Ey(i,j,k+2*dk) +  8*Ey(i,j,k+3*dk) -  Ey(i,j,k+4*dk) ) * odz60
        ad1_lE(2) = beta(1)*d1_f(1) + beta(2)*d1_f(2) + beta(3)*d1_f(3)

        d1_f(1) = di * (   2*Ez(i-2*di,j,k) - 24*Ez(i-di,j,k) - 35*Ez(i,j,k) + 80*Ez(i+di,j,k) &
                        - 30*Ez(i+2*di,j,k) +  8*Ez(i+3*di,j,k) -  Ez(i+4*di,j,k) ) * odx60
        d1_f(2) = dj * (   2*Ez(i,j-2*dj,k) - 24*Ez(i,j-dj,k) - 35*Ez(i,j,k) + 80*Ez(i,j+dj,k) &
                        - 30*Ez(i,j+2*dj,k) +  8*Ez(i,j+3*dj,k) -  Ez(i,j+4*dj,k) ) * ody60
        d1_f(3) = dk * (   2*Ez(i,j,k-2*dk) - 24*Ez(i,j,k-dk) - 35*Ez(i,j,k) + 80*Ez(i,j,k+dk) &
                        - 30*Ez(i,j,k+2*dk) +  8*Ez(i,j,k+3*dk) -  Ez(i,j,k+4*dk) ) * odz60
        ad1_lE(3) = beta(1)*d1_f(1) + beta(2)*d1_f(2) + beta(3)*d1_f(3)

        
        d1_f(1) = di * (   2*Ax(i-2*di,j,k) - 24*Ax(i-di,j,k) - 35*Ax(i,j,k) + 80*Ax(i+di,j,k) &
                        - 30*Ax(i+2*di,j,k) + 8*Ax(i+3*di,j,k) - Ax(i+4*di,j,k) ) * odx60
        d1_f(2) = dj * (   2*Ax(i,j-2*dj,k) - 24*Ax(i,j-dj,k) - 35*Ax(i,j,k) + 80*Ax(i,j+dj,k) &
                        - 30*Ax(i,j+2*dj,k) + 8*Ax(i,j+3*dj,k) - Ax(i,j+4*dj,k) ) * ody60
        d1_f(3) = dk * (   2*Ax(i,j,k-2*dk) - 24*Ax(i,j,k-dk) - 35*Ax(i,j,k) + 80*Ax(i,j,k+dk) &
                        - 30*Ax(i,j,k+2*dk) + 8*Ax(i,j,k+3*dk) - Ax(i,j,k+4*dk) ) * odz60
        ad1_lA(1) = beta(1)*d1_f(1) + beta(2)*d1_f(2) + beta(3)*d1_f(3)

        d1_f(1) = di * (   2*Ay(i-2*di,j,k) - 24*Ay(i-di,j,k) - 35*Ay(i,j,k) + 80*Ay(i+di,j,k) &
                        - 30*Ay(i+2*di,j,k) +  8*Ay(i+3*di,j,k) -  Ay(i+4*di,j,k) ) * odx60
        d1_f(2) = dj * (   2*Ay(i,j-2*dj,k) - 24*Ay(i,j-dj,k) - 35*Ay(i,j,k) + 80*Ay(i,j+dj,k) &
                        - 30*Ay(i,j+2*dj,k) +  8*Ay(i,j+3*dj,k) -  Ay(i,j+4*dj,k) ) * ody60
        d1_f(3) = dk * (   2*Ay(i,j,k-2*dk) - 24*Ay(i,j,k-dk) - 35*Ay(i,j,k) + 80*Ay(i,j,k+dk) &
                        - 30*Ay(i,j,k+2*dk) +  8*Ay(i,j,k+3*dk) -  Ay(i,j,k+4*dk) ) * odz60
        ad1_lA(2) = beta(1)*d1_f(1) + beta(2)*d1_f(2) + beta(3)*d1_f(3)

        d1_f(1) = di * (   2*Az(i-2*di,j,k) - 24*Az(i-di,j,k) - 35*Az(i,j,k) + 80*Az(i+di,j,k) &
                        - 30*Az(i+2*di,j,k) +  8*Az(i+3*di,j,k) -  Az(i+4*di,j,k) ) * odx60
        d1_f(2) = dj * (   2*Az(i,j-2*dj,k) - 24*Az(i,j-dj,k) - 35*Az(i,j,k) + 80*Az(i,j+dj,k) &
                        - 30*Az(i,j+2*dj,k) +  8*Az(i,j+3*dj,k) -  Az(i,j+4*dj,k) ) * ody60
        d1_f(3) = dk * (   2*Az(i,j,k-2*dk) - 24*Az(i,j,k-dk) - 35*Az(i,j,k) + 80*Az(i,j,k+dk) &
                        - 30*Az(i,j,k+2*dk) +  8*Az(i,j,k+3*dk) -  Az(i,j,k+4*dk) ) * odz60
        ad1_lA(3) = beta(1)*d1_f(1) + beta(2)*d1_f(2) + beta(3)*d1_f(3)

        
        d1_f(1) = di * (   2*Zeta(i-2*di,j,k) - 24*Zeta(i-di,j,k) - 35*Zeta(i,j,k) + 80*Zeta(i+di,j,k) &
                        - 30*Zeta(i+2*di,j,k) +  8*Zeta(i+3*di,j,k) -  Zeta(i+4*di,j,k) ) * odx60
        d1_f(2) = dj * (   2*Zeta(i,j-2*dj,k) - 24*Zeta(i,j-dj,k) - 35*Zeta(i,j,k) + 80*Zeta(i,j+dj,k) &
                        - 30*Zeta(i,j+2*dj,k) +  8*Zeta(i,j+3*dj,k) -  Zeta(i,j+4*dj,k) ) * ody60
        d1_f(3) = dk * (   2*Zeta(i,j,k-2*dk) - 24*Zeta(i,j,k-dk) - 35*Zeta(i,j,k) + 80*Zeta(i,j,k+dk) &
                        - 30*Zeta(i,j,k+2*dk) +  8*Zeta(i,j,k+3*dk) -  Zeta(i,j,k+4*dk) ) * odz60
        ad1_lZeta = beta(1)*d1_f(1) + beta(2)*d1_f(2) + beta(3)*d1_f(3)

        
        d1_f(1) = di * (   2*Aphi(i-2*di,j,k) - 24*Aphi(i-di,j,k) - 35*Aphi(i,j,k) + 80*Aphi(i+di,j,k) &
                        - 30*Aphi(i+2*di,j,k) +  8*Aphi(i+3*di,j,k) -  Aphi(i+4*di,j,k) ) * odx60
        d1_f(2) = dj * (   2*Aphi(i,j-2*dj,k) - 24*Aphi(i,j-dj,k) - 35*Aphi(i,j,k) + 80*Aphi(i,j+dj,k) &
                        - 30*Aphi(i,j+2*dj,k) +  8*Aphi(i,j+3*dj,k) -  Aphi(i,j+4*dj,k) ) * ody60
        d1_f(3) = dk * (   2*Aphi(i,j,k-2*dk) - 24*Aphi(i,j,k-dk) - 35*Aphi(i,j,k) + 80*Aphi(i,j,k+dk) &
                        - 30*Aphi(i,j,k+2*dk) +  8*Aphi(i,j,k+3*dk) -  Aphi(i,j,k+4*dk) ) * odz60
        ad1_lAphi = beta(1)*d1_f(1) + beta(2)*d1_f(2) + beta(3)*d1_f(3)

      else

        
        ad1_lE(1) = beta(1)*d1_lE(1,1) + beta(2)*d1_lE(1,2) + beta(3)*d1_lE(1,3)
        ad1_lE(2) = beta(1)*d1_lE(2,1) + beta(2)*d1_lE(2,2) + beta(3)*d1_lE(2,3)
        ad1_lE(3) = beta(1)*d1_lE(3,1) + beta(2)*d1_lE(3,2) + beta(3)*d1_lE(3,3)

        
        ad1_lA(1) = beta(1)*d1_lA(1,1) + beta(2)*d1_lA(1,2) + beta(3)*d1_lA(1,3)
        ad1_lA(2) = beta(1)*d1_lA(2,1) + beta(2)*d1_lA(2,2) + beta(3)*d1_lA(2,3)
        ad1_lA(3) = beta(1)*d1_lA(3,1) + beta(2)*d1_lA(3,2) + beta(3)*d1_lA(3,3)

        
        ad1_lZeta = beta(1)*d1_lZeta(1) + beta(2)*d1_lZeta(2) + beta(3)*d1_lZeta(3)

        
        ad1_lAphi = beta(1)*d1_lAphi(1) + beta(2)*d1_lAphi(2) + beta(3)*d1_lAphi(3)

      end if
      

    else
      call CCTK_Warn(0,812,"Proca_calc_rhs.F90","ProcaEvolve", "derivs_order not yet implemented.")
    end if

    
    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)
    


    
    cd_lA = d1_lA
    do a = 1, 3
      do b = 1, 3
        do m = 1, 3
          cd_lA(a,b) = cd_lA(a,b) - cf2(m,a,b) * lA(m)
       end do
      end do
    end do

    cd_dA = d2_lA
    do a = 1, 3
      do b = 1, 3
        do c = 1, 3
          do m = 1, 3
            cd_dA(a,b,c) = cd_dA(a,b,c) - cf2(m,a,c) * d1_lA(m,b)         &
                                        - cf2(m,b,c) * d1_lA(a,m)
           end do
         end do
       end do
     end do
    


    

    
    rhs_lE  = ad1_lE

    do a = 1, 3
      do m = 1, 3
        rhs_lE(a) = rhs_lE(a) - lE(m) * d1_beta(a,m)
      end do
    end do

    rhs_lE  = rhs_lE + alph * trk * lE

    do a = 1, 3
       do b = 1, 3
          rhs_lE(a) = rhs_lE(a) + alph * ch**conf_fac_exponent * hu(a,b) * d1_lZeta(b)      &
                                + alph * mu*mu * ch**conf_fac_exponent * hu(a,b) * lA(b)
          do c = 1, 3
             do m = 1, 3
                rhs_lE(a) = rhs_lE(a)                                                       &
                     + ch**(2*conf_fac_exponent) * hu(a,b) * hu(c,m) * d1_alph(m)           &
                                                 * ( d1_lA(c,b) - d1_lA(b,c) )              &
                     + alph * ch**(2*conf_fac_exponent) * hu(a,b) * hu(c,m)                 &
                                                 * ( cd_dA(m,b,c) - cd_dA(b,m,c) )          &
                     + 0.50d0*conf_fac_exponent * alph * ch**(2*conf_fac_exponent - 1)      &
                        * hu(a,b) * hu(c,m) * ( d1_lA(m,b)*d1_ch(c) - d1_lA(b,c)*d1_ch(m) )
             end do
          end do
       end do
    end do


    
    rhs_lA  = ad1_lA

    do a = 1, 3
      do m = 1, 3
        rhs_lA(a) = rhs_lA(a) + lA(m) * d1_beta(m,a)
      end do
    end do

    rhs_lA  = rhs_lA - alph * d1_lAphi - lAphi * d1_alph

    do a = 1,3
       do b = 1,3
          rhs_lA(a) = rhs_lA(a) - alph * hh(a,b) * lE(b) * ch**(-conf_fac_exponent)
       end do
    end do


    
    rhs_lAphi  = ad1_lAphi  + alph * trk * lAphi - alph * lZeta

    do a = 1, 3
        do b = 1, 3
           rhs_lAphi = rhs_lAphi                                                     &
                     + 0.5d0 * alph * conf_fac_exponent * ch**(conf_fac_exponent-1)  &
                             * hu(a,b) * lA(a) * d1_ch(b)                            &
                     - alph * ch**conf_fac_exponent * hu(a,b) * cd_lA(a,b)           &
                     - ch**conf_fac_exponent * hu(a,b) * lA(a) * d1_alph(b)
        end do
    end do


    
    rhs_lZeta = ad1_lZeta - alph * kappa * lZeta + alph * mu*mu * lAphi

    do a = 1, 3
       rhs_lZeta = rhs_lZeta + alph * d1_lE(a,a)                         &
                 - 1.5d0 * conf_fac_exponent * alph * lE(a) * d1_ch(a) / ch
    end do

    rhs_lZeta = Zeta_Omega_fac * rhs_lZeta

    

    
    
    
    
    
    
    
    
    
    


    

    rhs_Ex(i,j,k) = rhs_lE(1)
    rhs_Ey(i,j,k) = rhs_lE(2)
    rhs_Ez(i,j,k) = rhs_lE(3)

    rhs_Ax(i,j,k) = rhs_lA(1)
    rhs_Ay(i,j,k) = rhs_lA(2)
    rhs_Az(i,j,k) = rhs_lA(3)

    rhs_Zeta(i,j,k) = rhs_lZeta

    rhs_Aphi(i,j,k) = rhs_lAphi

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

end subroutine Proca_calc_rhs




subroutine Proca_calc_rhs_bdry( cctk_dim,cctk_gsh,cctk_lsh,cctk_lbnd,cctk_ubnd,cctk_tile_min,cctk_tile_max,cctk_ash,cctk_alignment,&
  &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_ash2,cctk&
  &_ash3, Aphi,Aphi_p,Aphi_p_p,Ax,Ax_p,Ax_p_p,Ay,Ay_p,Ay_p_p,Az,Az_p,Az_p_p,Ex,Ex_p,Ex_p_p,Ey,Ey_p,Ey_p_p,Ez,Ez_p,Ez_p_p,Zeta,Zeta_&
  &p,Zeta_p_p,alp,alp_p,alp_p_p,betax,betax_p,betax_p_p,betay,betay_p,betay_p_p,betaz,betaz_p,betaz_p_p,chi,coarse_dx,coarse_dy,coa&
  &rse_dz,conformal_state,dtalp,dtalp_p,dtalp_p_p,dtbetax,dtbetax_p,dtbetax_p_p,dtbetay,dtbetay_p,dtbetay_p_p,dtbetaz,dtbetaz_p,dtb&
  &etaz_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,eT&
  &xx_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,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,hxx,hxy,hxz,hyy,hyz,hzz,kxx,k&
  &xx_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,psi,psix,psixx,psixy,psix&
  &z,psiy,psiyy,psiyz,psiz,psizz,r,rhs_Aphi,rhs_Ax,rhs_Ay,rhs_Az,rhs_Ex,rhs_Ey,rhs_Ez,rhs_Zeta,shift_state,stress_energy_state,trac&
  &ek,x,y,z )

  implicit none

  INTEGER cctk_dim
  integer, parameter :: cctki_use_cctk_dim = kind(cctk_dim)
  INTEGER cctk_gsh (cctk_dim)
  integer, parameter :: cctki_use_cctk_gsh = kind(cctk_gsh)
  INTEGER cctk_lsh (cctk_dim)
  integer, parameter :: cctki_use_cctk_lsh = kind(cctk_lsh)
  INTEGER cctk_lbnd (cctk_dim)
  integer, parameter :: cctki_use_cctk_lbnd = kind(cctk_lbnd)
  INTEGER cctk_ubnd (cctk_dim)
  integer, parameter :: cctki_use_cctk_ubnd = kind(cctk_ubnd)
  INTEGER cctk_tile_min (cctk_dim)
  integer, parameter :: cctki_use_cctk_tile_min = kind(cctk_tile_min)
  INTEGER cctk_tile_max (cctk_dim)
  integer, parameter :: cctki_use_cctk_tile_max = kind(cctk_tile_max)
  INTEGER cctk_ash (cctk_dim)
  integer, parameter :: cctki_use_cctk_ash = kind(cctk_ash)
  INTEGER cctk_alignment
  integer, parameter :: cctki_use_cctk_alignment = kind(cctk_alignment)
  INTEGER cctk_alignment_offset
  integer, parameter :: cctki_use_cctk_alignment_offset = kind(cctk_alignment_offset)
  INTEGER cctk_from (cctk_dim)
  integer, parameter :: cctki_use_cctk_from = kind(cctk_from)
  INTEGER cctk_to (cctk_dim)
  integer, parameter :: cctki_use_cctk_to = kind(cctk_to)
  INTEGER cctk_bbox (2*cctk_dim)
  integer, parameter :: cctki_use_cctk_bbox = kind(cctk_bbox)
  REAL*8 cctk_delta_time
  integer, parameter :: cctki_use_cctk_delta_time = kind(cctk_delta_time)
  REAL*8 cctk_time
  integer, parameter :: cctki_use_cctk_time = kind(cctk_time)
  REAL*8 cctk_delta_space (cctk_dim)
  integer, parameter :: cctki_use_cctk_delta_space = kind(cctk_delta_space)
  REAL*8 cctk_origin_space (cctk_dim)
  integer, parameter :: cctki_use_cctk_origin_space = kind(cctk_origin_space)
  INTEGER cctk_levfac (cctk_dim)
  integer, parameter :: cctki_use_cctk_levfac = kind(cctk_levfac)
  INTEGER cctk_levoff (cctk_dim)
  integer, parameter :: cctki_use_cctk_levoff = kind(cctk_levoff)
  INTEGER cctk_levoffdenom (cctk_dim)
  integer, parameter :: cctki_use_cctk_levoffdenom = kind(cctk_levoffdenom)
  INTEGER cctk_timefac
  integer, parameter :: cctki_use_cctk_timefac = kind(cctk_timefac)
  INTEGER cctk_convlevel
  integer, parameter :: cctki_use_cctk_convlevel = kind(cctk_convlevel)
  INTEGER cctk_convfac
  integer, parameter :: cctki_use_cctk_convfac = kind(cctk_convfac)
  INTEGER cctk_nghostzones (cctk_dim)
  integer, parameter :: cctki_use_cctk_nghostzones = kind(cctk_nghostzones)
  INTEGER cctk_iteration
  integer, parameter :: cctki_use_cctk_iteration = kind(cctk_iteration)
  integer*8 cctkGH
  integer, parameter :: cctki_use_cctkGH = kind(cctkGH)
  INTEGER cctk_ash1
  integer, parameter :: cctki_use_cctk_ash1 = kind(cctk_ash1)
  INTEGER cctk_ash2
  integer, parameter :: cctki_use_cctk_ash2 = kind(cctk_ash2)
  INTEGER cctk_ash3
  integer, parameter :: cctki_use_cctk_ash3 = kind(cctk_ash3)
  REAL*8 Aphi (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Aphi = kind(Aphi)
  REAL*8 Aphi_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Aphi_p = kind(Aphi_p)
  REAL*8 Aphi_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Aphi_p_p = kind(Aphi_p_p)
  REAL*8 Ax (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Ax = kind(Ax)
  REAL*8 Ax_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Ax_p = kind(Ax_p)
  REAL*8 Ax_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Ax_p_p = kind(Ax_p_p)
  REAL*8 Ay (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Ay = kind(Ay)
  REAL*8 Ay_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Ay_p = kind(Ay_p)
  REAL*8 Ay_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Ay_p_p = kind(Ay_p_p)
  REAL*8 Az (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Az = kind(Az)
  REAL*8 Az_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Az_p = kind(Az_p)
  REAL*8 Az_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Az_p_p = kind(Az_p_p)
  REAL*8 Ex (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Ex = kind(Ex)
  REAL*8 Ex_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Ex_p = kind(Ex_p)
  REAL*8 Ex_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Ex_p_p = kind(Ex_p_p)
  REAL*8 Ey (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Ey = kind(Ey)
  REAL*8 Ey_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Ey_p = kind(Ey_p)
  REAL*8 Ey_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Ey_p_p = kind(Ey_p_p)
  REAL*8 Ez (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Ez = kind(Ez)
  REAL*8 Ez_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Ez_p = kind(Ez_p)
  REAL*8 Ez_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Ez_p_p = kind(Ez_p_p)
  REAL*8 Zeta (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Zeta = kind(Zeta)
  REAL*8 Zeta_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Zeta_p = kind(Zeta_p)
  REAL*8 Zeta_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_Zeta_p_p = kind(Zeta_p_p)
  REAL*8 alp (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_alp = kind(alp)
  REAL*8 alp_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_alp_p = kind(alp_p)
  REAL*8 alp_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_alp_p_p = kind(alp_p_p)
  REAL*8 betax (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_betax = kind(betax)
  REAL*8 betax_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_betax_p = kind(betax_p)
  REAL*8 betax_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_betax_p_p = kind(betax_p_p)
  REAL*8 betay (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_betay = kind(betay)
  REAL*8 betay_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_betay_p = kind(betay_p)
  REAL*8 betay_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_betay_p_p = kind(betay_p_p)
  REAL*8 betaz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_betaz = kind(betaz)
  REAL*8 betaz_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_betaz_p = kind(betaz_p)
  REAL*8 betaz_p_p (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_betaz_p_p = kind(betaz_p_p)
  REAL*8 chi (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_chi = kind(chi)
  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)
  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 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 hxx (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_hxx = kind(hxx)
  REAL*8 hxy (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_hxy = kind(hxy)
  REAL*8 hxz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_hxz = kind(hxz)
  REAL*8 hyy (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_hyy = kind(hyy)
  REAL*8 hyz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_hyz = kind(hyz)
  REAL*8 hzz (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_hzz = kind(hzz)
  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 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_Aphi (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_Aphi = kind(rhs_Aphi)
  REAL*8 rhs_Ax (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_Ax = kind(rhs_Ax)
  REAL*8 rhs_Ay (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_Ay = kind(rhs_Ay)
  REAL*8 rhs_Az (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_Az = kind(rhs_Az)
  REAL*8 rhs_Ex (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_Ex = kind(rhs_Ex)
  REAL*8 rhs_Ey (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_Ey = kind(rhs_Ey)
  REAL*8 rhs_Ez (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_Ez = kind(rhs_Ez)
  REAL*8 rhs_Zeta (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_rhs_Zeta = kind(rhs_Zeta)
  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 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 ProcaEvolve_MaxNumConstrainedVars
  integer, parameter :: cctki_use_ProcaEvolve_MaxNumConstrainedVars = kind(ProcaEvolve_MaxNumConstrainedVars)
  INTEGER*4 ProcaEvolve_MaxNumEvolvedVars
  integer, parameter :: cctki_use_ProcaEvolve_MaxNumEvolvedVars = kind(ProcaEvolve_MaxNumEvolvedVars)
  INTEGER*4 ProcaEvolve_MaxNumSandRVars
  integer, parameter :: cctki_use_ProcaEvolve_MaxNumSandRVars = kind(ProcaEvolve_MaxNumSandRVars)
  COMMON /ProcaEvolverest/ProcaEvolve_MaxNumConstrainedVars, ProcaEvolve_MaxNumEvolvedVars, ProcaEvolve_MaxNumSandRVars
  REAL*8 Zeta_Omega_fac
  integer, parameter :: cctki_use_Zeta_Omega_fac = kind(Zeta_Omega_fac)
  REAL*8 chi_floor
  integer, parameter :: cctki_use_chi_floor = kind(chi_floor)
  REAL*8 conf_fac_exponent
  integer, parameter :: cctki_use_conf_fac_exponent = kind(conf_fac_exponent)
  REAL*8 kappa
  integer, parameter :: cctki_use_kappa = kind(kappa)
  INTEGER*4 derivs_order
  integer, parameter :: cctki_use_derivs_order = kind(derivs_order)
  INTEGER*4 use_advection_stencils
  integer, parameter :: cctki_use_use_advection_stencils = kind(use_advection_stencils)
  COMMON /ProcaEvolvepriv/Zeta_Omega_fac, chi_floor, conf_fac_exponent, kappa, derivs_order, use_advection_stencils
  INTEGER*4 CCTKH0
  integer, parameter :: cctki_use_CCTKH0 = kind(CCTKH0)
  INTEGER*4 CCTKH1
  integer, parameter :: cctki_use_CCTKH1 = kind(CCTKH1)
  INTEGER*4 CCTKH2
  integer, parameter :: cctki_use_CCTKH2 = kind(CCTKH2)
  INTEGER*4 CCTKH3
  integer, parameter :: cctki_use_CCTKH3 = kind(CCTKH3)
  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 CCTKH4
  integer, parameter :: cctki_use_CCTKH4 = kind(CCTKH4)
  INTEGER*4 MoL_Num_SaveAndRestore_Vars
  integer, parameter :: cctki_use_MoL_Num_SaveAndRestore_Vars = kind(MoL_Num_SaveAndRestore_Vars)
  INTEGER*4 CCTKH5
  integer, parameter :: cctki_use_CCTKH5 = kind(CCTKH5)
  COMMON /METHODOFLINESrest/CCTKH0, CCTKH1, CCTKH2, CCTKH3, MoL_Num_Constrained_Vars, MoL_Num_Evolved_Vars, CCTKH4, MoL_Num_SaveAnd&
  &Restore_Vars, CCTKH5
  REAL*8 A0 (3)
  integer, parameter :: cctki_use_A0 = kind(A0)
  REAL*8 Aphi0
  integer, parameter :: cctki_use_Aphi0 = kind(Aphi0)
  REAL*8 E0 (3)
  integer, parameter :: cctki_use_E0 = kind(E0)
  REAL*8 mu
  integer, parameter :: cctki_use_mu = kind(mu)
  integer*8 evolution_method
  integer, parameter :: cctki_use_evolution_method = kind(evolution_method)
  integer*8 CCTKH6
  integer, parameter :: cctki_use_CCTKH6 = kind(CCTKH6)
  INTEGER*4 n_A (3)
  integer, parameter :: cctki_use_n_A = kind(n_A)
  INTEGER*4 n_Aphi
  integer, parameter :: cctki_use_n_Aphi = kind(n_Aphi)
  INTEGER*4 n_E (3)
  integer, parameter :: cctki_use_n_E = kind(n_E)
  INTEGER*4 n_Zeta
  integer, parameter :: cctki_use_n_Zeta = kind(n_Zeta)
  COMMON /PROCABASErest/A0, Aphi0, E0, mu, evolution_method, CCTKH6, n_A, n_Aphi, n_E, n_Zeta
  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 Boundary_SelectVarForBC (GH, faces, boundary_width, table_handle, var_name, bc_name)
  implicit none
  integer*8 GH
  INTEGER*4 faces
  INTEGER*4 boundary_width
  INTEGER*4 table_handle
  character(*) var_name
  character(*) bc_name
  end function Boundary_SelectVarForBC
  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 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, parameter :: one  = 1.0
  REAL*8, parameter :: zero = 0.0
  INTEGER*4 ierr

  ierr = NewRad_Apply(cctkGH, Ex, rhs_Ex, E0(1), one, n_E(1))
  ierr = NewRad_Apply(cctkGH, Ey, rhs_Ey, E0(2), one, n_E(2))
  ierr = NewRad_Apply(cctkGH, Ez, rhs_Ez, E0(3), one, n_E(3))

  ierr = NewRad_Apply(cctkGH, Ax, rhs_Ax, A0(1), one, n_A(1))
  ierr = NewRad_Apply(cctkGH, Ay, rhs_Ay, A0(2), one, n_A(2))
  ierr = NewRad_Apply(cctkGH, Az, rhs_Az, A0(3), one, n_A(3))

  ierr = NewRad_Apply(cctkGH, Aphi, rhs_Aphi, Aphi0, one, n_Aphi)
  ierr = NewRad_Apply(cctkGH, Zeta, rhs_Zeta, zero, one, n_Zeta)

end subroutine Proca_calc_rhs_bdry
