




















































































































 


































































































































































 

















































































































































































































































































































































 

























































































 










 








































































































































































































































 










                                
                                

                                
                                

                                
                                

                                




















































































subroutine Coordinates_SetInverseJacobian (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_l&
  &evfac,cctk_levoff,cctk_levoffdenom,cctk_timefac,cctk_convlevel,cctk_convfac,cctk_nghostzones,cctk_iteration,cctkGH,cctk_ash1,cct&
  &k_ash2,cctk_ash3, J11,J12,J13,J21,J22,J23,J31,J32,J33,coarse_dx,coarse_dy,coarse_dz,dJ111,dJ112,dJ113,dJ122,dJ123,dJ133,dJ211,dJ&
  &212,dJ213,dJ222,dJ223,dJ233,dJ311,dJ312,dJ313,dJ322,dJ323,dJ333,general_coordinates,iJ11,iJ12,iJ13,iJ21,iJ22,iJ23,iJ31,iJ32,iJ33&
  &,interpolate_boundary_points,inverse_jacobian_state,jacobian_derivative_state,jacobian_state,r,volume_form,volume_form_state,x,y&
  &,z)

  use matdet
  use matinv
  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 J11 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_J11 = kind(J11)
  REAL*8 J12 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_J12 = kind(J12)
  REAL*8 J13 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_J13 = kind(J13)
  REAL*8 J21 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_J21 = kind(J21)
  REAL*8 J22 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_J22 = kind(J22)
  REAL*8 J23 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_J23 = kind(J23)
  REAL*8 J31 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_J31 = kind(J31)
  REAL*8 J32 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_J32 = kind(J32)
  REAL*8 J33 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_J33 = kind(J33)
  REAL*8 coarse_dx
  integer, parameter :: cctki_use_coarse_dx = kind(coarse_dx)
  REAL*8 coarse_dy
  integer, parameter :: cctki_use_coarse_dy = kind(coarse_dy)
  REAL*8 coarse_dz
  integer, parameter :: cctki_use_coarse_dz = kind(coarse_dz)
  REAL*8 dJ111 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_dJ111 = kind(dJ111)
  REAL*8 dJ112 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_dJ112 = kind(dJ112)
  REAL*8 dJ113 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_dJ113 = kind(dJ113)
  REAL*8 dJ122 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_dJ122 = kind(dJ122)
  REAL*8 dJ123 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_dJ123 = kind(dJ123)
  REAL*8 dJ133 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_dJ133 = kind(dJ133)
  REAL*8 dJ211 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_dJ211 = kind(dJ211)
  REAL*8 dJ212 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_dJ212 = kind(dJ212)
  REAL*8 dJ213 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_dJ213 = kind(dJ213)
  REAL*8 dJ222 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_dJ222 = kind(dJ222)
  REAL*8 dJ223 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_dJ223 = kind(dJ223)
  REAL*8 dJ233 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_dJ233 = kind(dJ233)
  REAL*8 dJ311 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_dJ311 = kind(dJ311)
  REAL*8 dJ312 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_dJ312 = kind(dJ312)
  REAL*8 dJ313 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_dJ313 = kind(dJ313)
  REAL*8 dJ322 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_dJ322 = kind(dJ322)
  REAL*8 dJ323 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_dJ323 = kind(dJ323)
  REAL*8 dJ333 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_dJ333 = kind(dJ333)
  INTEGER*4 general_coordinates
  integer, parameter :: cctki_use_general_coordinates = kind(general_coordinates)
  REAL*8 iJ11 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_iJ11 = kind(iJ11)
  REAL*8 iJ12 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_iJ12 = kind(iJ12)
  REAL*8 iJ13 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_iJ13 = kind(iJ13)
  REAL*8 iJ21 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_iJ21 = kind(iJ21)
  REAL*8 iJ22 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_iJ22 = kind(iJ22)
  REAL*8 iJ23 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_iJ23 = kind(iJ23)
  REAL*8 iJ31 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_iJ31 = kind(iJ31)
  REAL*8 iJ32 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_iJ32 = kind(iJ32)
  REAL*8 iJ33 (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_iJ33 = kind(iJ33)
  INTEGER*4 interpolate_boundary_points
  integer, parameter :: cctki_use_interpolate_boundary_points = kind(interpolate_boundary_points)
  INTEGER*4 inverse_jacobian_state
  integer, parameter :: cctki_use_inverse_jacobian_state = kind(inverse_jacobian_state)
  INTEGER*4 jacobian_derivative_state
  integer, parameter :: cctki_use_jacobian_derivative_state = kind(jacobian_derivative_state)
  INTEGER*4 jacobian_state
  integer, parameter :: cctki_use_jacobian_state = kind(jacobian_state)
  REAL*8 r (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_r = kind(r)
  REAL*8 volume_form (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_volume_form = kind(volume_form)
  INTEGER*4 volume_form_state
  integer, parameter :: cctki_use_volume_form_state = kind(volume_form_state)
  REAL*8 x (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_x = kind(x)
  REAL*8 y (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_y = kind(y)
  REAL*8 z (cctk_ash1,cctk_ash2,cctk_ash3)
  integer, parameter :: cctki_use_z = kind(z)
  integer*8 coordinate_system
  integer, parameter :: cctki_use_coordinate_system = kind(coordinate_system)
  COMMON /Coordinatesrest/coordinate_system
  REAL*8 box_radius
  integer, parameter :: cctki_use_box_radius = kind(box_radius)
  REAL*8 cylinder_inner_radius
  integer, parameter :: cctki_use_cylinder_inner_radius = kind(cylinder_inner_radius)
  REAL*8 cylinder_outer_radius
  integer, parameter :: cctki_use_cylinder_outer_radius = kind(cylinder_outer_radius)
  REAL*8 cylinder_radius
  integer, parameter :: cctki_use_cylinder_radius = kind(cylinder_radius)
  REAL*8 cylinder_zmax
  integer, parameter :: cctki_use_cylinder_zmax = kind(cylinder_zmax)
  REAL*8 cylinder_zmin
  integer, parameter :: cctki_use_cylinder_zmin = kind(cylinder_zmin)
  REAL*8 h_cartesian
  integer, parameter :: cctki_use_h_cartesian = kind(h_cartesian)
  REAL*8 h_radial
  integer, parameter :: cctki_use_h_radial = kind(h_radial)
  REAL*8 h_radial_1
  integer, parameter :: cctki_use_h_radial_1 = kind(h_radial_1)
  REAL*8 h_radial_inner
  integer, parameter :: cctki_use_h_radial_inner = kind(h_radial_inner)
  REAL*8 h_radial_outer
  integer, parameter :: cctki_use_h_radial_outer = kind(h_radial_outer)
  REAL*8 h_z
  integer, parameter :: cctki_use_h_z = kind(h_z)
  REAL*8 patch_one_xmax
  integer, parameter :: cctki_use_patch_one_xmax = kind(patch_one_xmax)
  REAL*8 patch_one_xmin
  integer, parameter :: cctki_use_patch_one_xmin = kind(patch_one_xmin)
  REAL*8 patch_one_ymax
  integer, parameter :: cctki_use_patch_one_ymax = kind(patch_one_ymax)
  REAL*8 patch_one_ymin
  integer, parameter :: cctki_use_patch_one_ymin = kind(patch_one_ymin)
  REAL*8 patch_one_zmax
  integer, parameter :: cctki_use_patch_one_zmax = kind(patch_one_zmax)
  REAL*8 patch_one_zmin
  integer, parameter :: cctki_use_patch_one_zmin = kind(patch_one_zmin)
  REAL*8 patch_two_xmax
  integer, parameter :: cctki_use_patch_two_xmax = kind(patch_two_xmax)
  REAL*8 patch_two_xmin
  integer, parameter :: cctki_use_patch_two_xmin = kind(patch_two_xmin)
  REAL*8 patch_two_ymax
  integer, parameter :: cctki_use_patch_two_ymax = kind(patch_two_ymax)
  REAL*8 patch_two_ymin
  integer, parameter :: cctki_use_patch_two_ymin = kind(patch_two_ymin)
  REAL*8 patch_two_zmax
  integer, parameter :: cctki_use_patch_two_zmax = kind(patch_two_zmax)
  REAL*8 patch_two_zmin
  integer, parameter :: cctki_use_patch_two_zmin = kind(patch_two_zmin)
  REAL*8 patch_xmax
  integer, parameter :: cctki_use_patch_xmax = kind(patch_xmax)
  REAL*8 patch_xmin
  integer, parameter :: cctki_use_patch_xmin = kind(patch_xmin)
  REAL*8 patch_ymax
  integer, parameter :: cctki_use_patch_ymax = kind(patch_ymax)
  REAL*8 patch_ymin
  integer, parameter :: cctki_use_patch_ymin = kind(patch_ymin)
  REAL*8 patch_zmax
  integer, parameter :: cctki_use_patch_zmax = kind(patch_zmax)
  REAL*8 patch_zmin
  integer, parameter :: cctki_use_patch_zmin = kind(patch_zmin)
  REAL*8 sphere_inner_radius
  integer, parameter :: cctki_use_sphere_inner_radius = kind(sphere_inner_radius)
  REAL*8 sphere_medium_radius
  integer, parameter :: cctki_use_sphere_medium_radius = kind(sphere_medium_radius)
  REAL*8 sphere_outer_radius
  integer, parameter :: cctki_use_sphere_outer_radius = kind(sphere_outer_radius)
  REAL*8 stretch_rmax_1
  integer, parameter :: cctki_use_stretch_rmax_1 = kind(stretch_rmax_1)
  REAL*8 stretch_rmin_1
  integer, parameter :: cctki_use_stretch_rmin_1 = kind(stretch_rmin_1)
  REAL*8 theta_min
  integer, parameter :: cctki_use_theta_min = kind(theta_min)
  REAL*8 transition_radius
  integer, parameter :: cctki_use_transition_radius = kind(transition_radius)
  integer*8 symmetry
  integer, parameter :: cctki_use_symmetry = kind(symmetry)
  INTEGER*4 MonteCarloSeed
  integer, parameter :: cctki_use_MonteCarloSeed = kind(MonteCarloSeed)
  INTEGER*4 additional_overlap_size
  integer, parameter :: cctki_use_additional_overlap_size = kind(additional_overlap_size)
  INTEGER*4 additional_symmetry_size
  integer, parameter :: cctki_use_additional_symmetry_size = kind(additional_symmetry_size)
  INTEGER*4 cubical_inner_boundary
  integer, parameter :: cctki_use_cubical_inner_boundary = kind(cubical_inner_boundary)
  INTEGER*4 internal_outer_boundaries
  integer, parameter :: cctki_use_internal_outer_boundaries = kind(internal_outer_boundaries)
  INTEGER*4 nMonteCarloParticles
  integer, parameter :: cctki_use_nMonteCarloParticles = kind(nMonteCarloParticles)
  INTEGER*4 n_angular
  integer, parameter :: cctki_use_n_angular = kind(n_angular)
  INTEGER*4 n_angular_inner
  integer, parameter :: cctki_use_n_angular_inner = kind(n_angular_inner)
  INTEGER*4 n_angular_outer
  integer, parameter :: cctki_use_n_angular_outer = kind(n_angular_outer)
  INTEGER*4 n_angular_phi
  integer, parameter :: cctki_use_n_angular_phi = kind(n_angular_phi)
  INTEGER*4 n_angular_theta
  integer, parameter :: cctki_use_n_angular_theta = kind(n_angular_theta)
  INTEGER*4 n_xy
  integer, parameter :: cctki_use_n_xy = kind(n_xy)
  INTEGER*4 ncells_x
  integer, parameter :: cctki_use_ncells_x = kind(ncells_x)
  INTEGER*4 ncells_y
  integer, parameter :: cctki_use_ncells_y = kind(ncells_y)
  INTEGER*4 ncells_z
  integer, parameter :: cctki_use_ncells_z = kind(ncells_z)
  INTEGER*4 outer_boundary_size
  integer, parameter :: cctki_use_outer_boundary_size = kind(outer_boundary_size)
  INTEGER*4 patch_boundary_size
  integer, parameter :: cctki_use_patch_boundary_size = kind(patch_boundary_size)
  INTEGER*4 patch_one_ncells_x
  integer, parameter :: cctki_use_patch_one_ncells_x = kind(patch_one_ncells_x)
  INTEGER*4 patch_one_ncells_y
  integer, parameter :: cctki_use_patch_one_ncells_y = kind(patch_one_ncells_y)
  INTEGER*4 patch_one_ncells_z
  integer, parameter :: cctki_use_patch_one_ncells_z = kind(patch_one_ncells_z)
  INTEGER*4 patch_two_ncells_x
  integer, parameter :: cctki_use_patch_two_ncells_x = kind(patch_two_ncells_x)
  INTEGER*4 patch_two_ncells_y
  integer, parameter :: cctki_use_patch_two_ncells_y = kind(patch_two_ncells_y)
  INTEGER*4 patch_two_ncells_z
  integer, parameter :: cctki_use_patch_two_ncells_z = kind(patch_two_ncells_z)
  INTEGER*4 radial_stretch
  integer, parameter :: cctki_use_radial_stretch = kind(radial_stretch)
  INTEGER*4 register_symmetry
  integer, parameter :: cctki_use_register_symmetry = kind(register_symmetry)
  INTEGER*4 shiftout_outer_boundaries
  integer, parameter :: cctki_use_shiftout_outer_boundaries = kind(shiftout_outer_boundaries)
  INTEGER*4 stagger_outer_boundaries
  integer, parameter :: cctki_use_stagger_outer_boundaries = kind(stagger_outer_boundaries)
  INTEGER*4 stagger_patch_boundaries
  integer, parameter :: cctki_use_stagger_patch_boundaries = kind(stagger_patch_boundaries)
  INTEGER*4 store_inverse_jacobian
  integer, parameter :: cctki_use_store_inverse_jacobian = kind(store_inverse_jacobian)
  INTEGER*4 store_jacobian
  integer, parameter :: cctki_use_store_jacobian = kind(store_jacobian)
  INTEGER*4 store_jacobian_derivative
  integer, parameter :: cctki_use_store_jacobian_derivative = kind(store_jacobian_derivative)
  INTEGER*4 store_volume_form
  integer, parameter :: cctki_use_store_volume_form = kind(store_volume_form)
  INTEGER*4 verbose
  integer, parameter :: cctki_use_verbose = kind(verbose)
  COMMON /Coordinatespriv/box_radius, cylinder_inner_radius, cylinder_outer_radius, cylinder_radius, cylinder_zmax, cylinder_zmin, &
  &h_cartesian, h_radial, h_radial_1, h_radial_inner, h_radial_outer, h_z, patch_one_xmax, patch_one_xmin, patch_one_ymax, patch_on&
  &e_ymin, patch_one_zmax, patch_one_zmin, patch_two_xmax, patch_two_xmin, patch_two_ymax, patch_two_ymin, patch_two_zmax, patch_tw&
  &o_zmin, patch_xmax, patch_xmin, patch_ymax, patch_ymin, patch_zmax, patch_zmin, sphere_inner_radius, sphere_medium_radius, spher&
  &e_outer_radius, stretch_rmax_1, stretch_rmin_1, theta_min, transition_radius, symmetry, MonteCarloSeed, additional_overlap_size,&
  & additional_symmetry_size, cubical_inner_boundary, internal_outer_boundaries, nMonteCarloParticles, n_angular, n_angular_inner, &
  &n_angular_outer, n_angular_phi, n_angular_theta, n_xy, ncells_x, ncells_y, ncells_z, outer_boundary_size, patch_boundary_size, p&
  &atch_one_ncells_x, patch_one_ncells_y, patch_one_ncells_z, patch_two_ncells_x, patch_two_ncells_y, patch_two_ncells_z, radial_st&
  &retch, register_symmetry, shiftout_outer_boundaries, stagger_outer_boundaries, stagger_patch_boundaries, store_inverse_jacobian,&
  & store_jacobian, store_jacobian_derivative, store_volume_form, verbose
  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 CoordinatesSymmetry_Interpolate (cctkGH, N_dims, local_interp_handle, param_table_handle, coord_system_handle,&
  & N_interp_points, interp_coords_type, interp_coords, N_input_arrays, input_array_indices, N_output_arrays, output_array_types, o&
  &utput_arrays, faces)
  implicit none
  integer*8 cctkGH
  INTEGER*4 N_dims
  INTEGER*4 local_interp_handle
  INTEGER*4 param_table_handle
  INTEGER*4 coord_system_handle
  INTEGER*4 N_interp_points
  INTEGER*4 interp_coords_type
  integer*8 interp_coords(*)
  INTEGER*4 N_input_arrays
  INTEGER*4 input_array_indices(*)
  INTEGER*4 N_output_arrays
  INTEGER*4 output_array_types(*)
  integer*8 output_arrays(*)
  INTEGER*4 faces
  end function CoordinatesSymmetry_Interpolate
  end interface
  interface
  INTEGER*4 function GetBoundarySpecification (size, nboundaryzones, is_internal, is_staggered, shiftout)
  implicit none
  INTEGER*4 size
  INTEGER*4 nboundaryzones(*)
  INTEGER*4 is_internal(*)
  INTEGER*4 is_staggered(*)
  INTEGER*4 shiftout(*)
  end function GetBoundarySpecification
  end interface
  interface
  INTEGER*4 function GetRefinementLevel (cctkGH)
  implicit none
  integer*8 cctkGH
  end function GetRefinementLevel
  end interface
  interface
  INTEGER*4 function MultiPatch_ConvertFromPhysicalBoundary (map, size, physical_min, physical_max, interior_min, interior_max, ext&
  &erior_min, exterior_max, spacing)
  implicit none
  INTEGER*4 map
  INTEGER*4 size
  REAL*8 physical_min(*)
  REAL*8 physical_max(*)
  REAL*8 interior_min(*)
  REAL*8 interior_max(*)
  REAL*8 exterior_min(*)
  REAL*8 exterior_max(*)
  REAL*8 spacing(*)
  end function MultiPatch_ConvertFromPhysicalBoundary
  end interface
  interface
  INTEGER*4 function MultiPatch_GetBbox (cctkGH, size, bbox)
  implicit none
  integer*8 cctkGH
  INTEGER*4 size
  INTEGER*4 bbox(*)
  end function MultiPatch_GetBbox
  end interface
  interface
  INTEGER*4 function MultiPatch_GetBoundarySpecification (map, size, nboundaryzones, is_internal, is_staggered, shiftout)
  implicit none
  INTEGER*4 map
  INTEGER*4 size
  INTEGER*4 nboundaryzones(*)
  INTEGER*4 is_internal(*)
  INTEGER*4 is_staggered(*)
  INTEGER*4 shiftout(*)
  end function MultiPatch_GetBoundarySpecification
  end interface
  interface
  INTEGER*4 function MultiPatch_GetDomainSpecification (map, size, physical_min, physical_max, interior_min, interior_max, exterior&
  &_min, exterior_max, spacing)
  implicit none
  INTEGER*4 map
  INTEGER*4 size
  REAL*8 physical_min(*)
  REAL*8 physical_max(*)
  REAL*8 interior_min(*)
  REAL*8 interior_max(*)
  REAL*8 exterior_min(*)
  REAL*8 exterior_max(*)
  REAL*8 spacing(*)
  end function MultiPatch_GetDomainSpecification
  end interface
  interface
  INTEGER*4 function MultiPatch_GetInnerRadius ()
  implicit none
  end function MultiPatch_GetInnerRadius
  end interface
  interface
  INTEGER*4 function MultiPatch_GetMap (cctkGH)
  implicit none
  integer*8 cctkGH
  end function MultiPatch_GetMap
  end interface
  interface
  INTEGER*4 function MultiPatch_GetMaps (cctkGH)
  implicit none
  integer*8 cctkGH
  end function MultiPatch_GetMaps
  end interface
  interface
  INTEGER*4 function MultiPatch_GetSymmetryBoundaries (cctkGH, size, symmbnd)
  implicit none
  integer*8 cctkGH
  INTEGER*4 size
  INTEGER*4 symmbnd(*)
  end function MultiPatch_GetSymmetryBoundaries
  end interface
  interface
  INTEGER*4 function MultiPatch_GetSymmetrySpecification (map, size, type, symmetry_overlap, symmetry_dir)
  implicit none
  INTEGER*4 map
  INTEGER*4 size
  INTEGER*4 type(*)
  INTEGER*4 symmetry_overlap
  INTEGER*4 symmetry_dir
  end function MultiPatch_GetSymmetrySpecification
  end interface
  interface
  INTEGER*4 function MultiPatch_GetSystemSpecification (maps)
  implicit none
  INTEGER*4 maps
  end function MultiPatch_GetSystemSpecification
  end interface
  interface
  INTEGER*4 function MultiPatch_GlobalToLocal (cctkGH, ndims, npoints, globalcoords, patch, localcoords, dadx, ddadxdx)
  implicit none
  integer*8 cctkGH
  INTEGER*4 ndims
  INTEGER*4 npoints
  integer*8 globalcoords
  INTEGER*4 patch(*)
  integer*8 localcoords
  integer*8 dadx
  integer*8 ddadxdx
  end function MultiPatch_GlobalToLocal
  end interface
  interface
  INTEGER*4 function MultiPatch_LocalToGlobal (cctkGH, ndims, npoints, patch, localcoords, globalcoords, dxda, det_dxda, dadx, ddxd&
  &ada, ddadxdx, dddxdadada)
  implicit none
  integer*8 cctkGH
  INTEGER*4 ndims
  INTEGER*4 npoints
  INTEGER*4 patch(*)
  integer*8 localcoords
  integer*8 globalcoords
  integer*8 dxda
  integer*8 det_dxda
  integer*8 dadx
  integer*8 ddxdada
  integer*8 ddadxdx
  integer*8 dddxdadada
  end function MultiPatch_LocalToGlobal
  end interface
  interface
  INTEGER*4 function MultiPatch_MapIsCartesian (map)
  implicit none
  INTEGER*4 map
  end function MultiPatch_MapIsCartesian
  end interface
  interface
  INTEGER*4 function MultiPatch_ProvidesThornburg04 ()
  implicit none
  end function MultiPatch_ProvidesThornburg04
  end interface
  interface
  INTEGER*4 function MultiPatch_SymmetryInterpolateFaces (cctkGH, N_dims, local_interp_handle, param_table_handle, coord_system_han&
  &dle, N_interp_points, interp_coords_type, interp_coords, N_input_arrays, input_array_indices, N_output_arrays, output_array_type&
  &s, output_arrays, faces)
  implicit none
  integer*8 cctkGH
  INTEGER*4 N_dims
  INTEGER*4 local_interp_handle
  INTEGER*4 param_table_handle
  INTEGER*4 coord_system_handle
  INTEGER*4 N_interp_points
  INTEGER*4 interp_coords_type
  integer*8 interp_coords(*)
  INTEGER*4 N_input_arrays
  INTEGER*4 input_array_indices(*)
  INTEGER*4 N_output_arrays
  INTEGER*4 output_array_types(*)
  integer*8 output_arrays(*)
  INTEGER*4 faces
  end function MultiPatch_SymmetryInterpolateFaces
  end interface
  interface
  INTEGER*4 function SymmetryInterpolateFaces (cctkGH, N_dims, local_interp_handle, param_table_handle, coord_system_handle, N_inte&
  &rp_points, interp_coords_type, interp_coords, N_input_arrays, input_array_indices, N_output_arrays, output_array_types, output_a&
  &rrays, faces)
  implicit none
  integer*8 cctkGH
  INTEGER*4 N_dims
  INTEGER*4 local_interp_handle
  INTEGER*4 param_table_handle
  INTEGER*4 coord_system_handle
  INTEGER*4 N_interp_points
  INTEGER*4 interp_coords_type
  integer*8 interp_coords(*)
  INTEGER*4 N_input_arrays
  INTEGER*4 input_array_indices(*)
  INTEGER*4 N_output_arrays
  INTEGER*4 output_array_types(*)
  integer*8 output_arrays(*)
  INTEGER*4 faces
  end function SymmetryInterpolateFaces
  end interface
  interface
  INTEGER*4 function SymmetryRegister (sym_name)
  implicit none
  character(*) sym_name
  end function SymmetryRegister
  end interface
  interface
  INTEGER*4 function SymmetryRegisterGrid (cctkGH, sym_handle, which_faces, symmetry_zone_width)
  implicit none
  integer*8 cctkGH
  INTEGER*4 sym_handle
  INTEGER*4 which_faces(*)
  INTEGER*4 symmetry_zone_width(*)
  end function SymmetryRegisterGrid
  end interface
  interface
  INTEGER*4 function SymmetryRegisterGridInterpolator (cctkGH, sym_handle, symmetry_interpolate)
  implicit none
  integer*8 cctkGH
  INTEGER*4 sym_handle
  external symmetry_interpolate
  INTEGER*4 symmetry_interpolate
  end function SymmetryRegisterGridInterpolator
  end interface
  interface
  INTEGER*4 function SymmetryTableHandleForGrid (cctkGH)
  implicit none
  integer*8 cctkGH
  end function SymmetryTableHandleForGrid
  end interface
  interface
  INTEGER*4 function Util_GetGridRanges (GH, istart, iend)
  implicit none
  integer*8 GH
  INTEGER*4 istart(*)
  INTEGER*4 iend(*)
  end function Util_GetGridRanges
  end interface
  interface
  subroutine Util_GetPhysicalBoundaries (GH, do_bc)
  implicit none
  integer*8 GH
  INTEGER*4 do_bc(*)
  end subroutine Util_GetPhysicalBoundaries
  end interface

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

  REAL*8                 :: detJ
  REAL*8, dimension(3,3) :: Jac, iJac

  volume_form_state      = store_volume_form
  inverse_Jacobian_state = store_inverse_jacobian

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

  !$omp parallel do private(i,j,k, detJ, Jac, iJac)
  do k = 1, nz
     do j = 1, ny
        do i = 1, nx

           Jac(1,1) = J11(i,j,k)
           Jac(1,2) = J12(i,j,k)
           Jac(1,3) = J13(i,j,k)
           Jac(2,1) = J21(i,j,k)
           Jac(2,2) = J22(i,j,k)
           Jac(2,3) = J23(i,j,k)
           Jac(3,1) = J31(i,j,k)
           Jac(3,2) = J32(i,j,k)
           Jac(3,3) = J33(i,j,k)

           if (store_volume_form /= 0) then
              call calc_det3 (Jac, detJ)

              
              
              
              volume_form(i,j,k) = detJ
           end if

           if (store_inverse_jacobian /= 0) then
              call calc_inv3 (Jac, iJac)

              iJ11(i,j,k) = iJac(1,1)
              iJ12(i,j,k) = iJac(1,2)
              iJ13(i,j,k) = iJac(1,3)
              iJ21(i,j,k) = iJac(2,1)
              iJ22(i,j,k) = iJac(2,2)
              iJ23(i,j,k) = iJac(2,3)
              iJ31(i,j,k) = iJac(3,1)
              iJ32(i,j,k) = iJac(3,2)
              iJ33(i,j,k) = iJac(3,3)
           end if

        end do
     end do
  end do

end subroutine Coordinates_SetInverseJacobian
