diff --git a/Jobs/NHSphere/JobNHSchaerHillDrySphere b/Jobs/NHSphere/JobNHSchaerHillDrySphere
index 774a519..484b712 100755
--- a/Jobs/NHSphere/JobNHSchaerHillDrySphere
+++ b/Jobs/NHSphere/JobNHSchaerHillDrySphere
@@ -1,4 +1,4 @@
-mpirun -n 6 julia --project Examples/testNHSphere.jl \
+mpirun -n 1 julia --project Examples/testNHSphere.jl \
   --Problem="SchaerSphericalSphere" \
   --Device="CPU" \
   --GPUType="Metal" \
diff --git a/src/Outputs/vtkOutputKernel.jl b/src/Outputs/vtkOutputKernel.jl
index f22fdbc..8988a0c 100644
--- a/src/Outputs/vtkOutputKernel.jl
+++ b/src/Outputs/vtkOutputKernel.jl
@@ -73,6 +73,27 @@ end
   end
 end
 
+@kernel function InterpolateCGDim2Kernel!(cCell,@Const(c),@Const(Inter),@Const(Glob),
+  ::Val{BANK}=Val(1)) where BANK
+  I, J   = @index(Local,  NTuple)
+  _,_,IF = @index(Global,  NTuple)
+
+
+  ColumnTilesDim = @uniform @groupsize()[3]
+  NF = @uniform @ndrange()[3]
+
+  @uniform N = size(Inter,3)
+
+  if IF <= NF
+    @inbounds cCell[I,J,IF] = eltype(cCell)(0)
+    for jP = 1 : N
+      for iP = 1 : N
+        @inbounds  cCell[I,J,IF] += Inter[I,J,iP,jP] * c[I,J,IF]
+      end
+    end
+  end
+end
+
 @kernel function InterpolateRhoKernel!(cCell,@Const(c),@Const(RhoC),@Const(Inter),@Const(Glob),
   ::Val{BANK}=Val(1)) where BANK
   I, J, iz   = @index(Local,  NTuple)
@@ -211,6 +232,25 @@ function InterpolateGPU!(cCell,c,Inter,Glob)
 
 end
 
+function InterpolateCGDim2GPU!(cCell,c,Inter,Glob)
+
+  backend = get_backend(c)
+  FT = eltype(c)
+
+  OrdPrint = size(Inter,1)
+  NF = size(Glob,2)
+
+# Ranges
+  NFG = min(div(256,OrdPrint*OrdPrint),NF)
+  group = (OrdPrint, OrdPrint, NFG)
+  ndrange = (OrdPrint, OrdPrint, NF)
+
+  KInterpolateCGDim2Kernel! = InterpolateCGDim2Kernel!(backend,group)
+  KInterpolateCGDim2Kernel!(cCell,c,Inter,Glob,ndrange=ndrange)
+  KernelAbstractions.synchronize(backend)
+
+end
+
 function InterpolateRhoGPU!(cCell,c,Rho,Inter,Glob)
 
   backend = get_backend(c)
diff --git a/src/Outputs/vtkSphere.jl b/src/Outputs/vtkSphere.jl
index 5fba6c8..ec00578 100644
--- a/src/Outputs/vtkSphere.jl
+++ b/src/Outputs/vtkSphere.jl
@@ -684,7 +684,12 @@ function unstructured_vtkOrography(Height,vtkGrid, NF, CG,  part::Int, nparts::I
   vtk_filename_noext = filename
   vtk = pvtk_grid(vtk_filename_noext, pts, cells; compress=3, part = part, nparts = nparts)
   HeightCell = zeros(OrdPrint*OrdPrint*NF) 
-  @views InterpolateCG!(HeightCell,Height,vtkInter,OrdPoly,OrdPrint,CG.Glob,NF,nz)
+  backend = get_backend(Height)
+  FTB = eltype(Height)
+  cCell = KernelAbstractions.zeros(backend,FTB,OrdPrint,OrdPrint,NF)
+  #@views InterpolateCG!(HeightCell,Height,vtkInter,OrdPoly,OrdPrint,CG.Glob,NF,nz)
+  InterpolateCGDim2GPU!(cCell,Height,vtkInter,CG.Glob)
+  copyto!(HeightCell,cCell)
   vtk["Height", VTKCellData()] = HeightCell
   outfiles=vtk_save(vtk);
   return outfiles::Vector{String}