diff --git a/Project.toml b/Project.toml
index 4d6eb29..68e77f8 100644
--- a/Project.toml
+++ b/Project.toml
@@ -26,7 +26,7 @@ WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 [compat]
 CUDA = "4 - 5"
 CairoMakie = "0.8 - 0.20"
-GeophysicalModelGenerator = "0.3 - 0.6"
+GeophysicalModelGenerator = "0.7"
 GeoParams = "0.5"
 Interpolations = "0.13 - 0.15"
 JLD2 = "0.4"
diff --git a/README.md b/README.md
index ea9c396..21c8d86 100644
--- a/README.md
+++ b/README.md
@@ -10,10 +10,14 @@ This easy to use and versatile package simulates the thermal evolution of magmat
 Below we give a number of example scripts that show how it can be used to simulate a number of scenarios.
 
 ## Contents
+- [MagmaThermoKinematics.jl](#magmathermokinematicsjl)
+  - [Contents](#contents)
   - [100-lines 2D example](#100-lines-2d-example)
   - [100-lines 3D example](#100-lines-3d-example)
   - [Dependencies](#dependencies)
   - [Installation](#installation)
+  - [Ongoing development](#ongoing-development)
+  - [Related work](#related-work)
 
 ## 100-lines 2D example
 A simple example that simulates the emplacement of dikes within the crust over a period of 10'000 years is shown below. 
@@ -41,7 +45,7 @@ Num                     =   Numeric_params(verbose=false)                   # No
 MatParam                =   (
         SetMaterialParams(Name="Rock", Phase=1, 
              Density    = ConstantDensity(ρ=2800kg/m^3),               
-           HeatCapacity = ConstantHeatCapacity(cp=1050J/kg/K),
+           HeatCapacity = ConstantHeatCapacity(Cp=1050J/kg/K),
            Conductivity = ConstantConductivity(k=1.5Watt/K/m),       
              LatentHeat = ConstantLatentHeat(Q_L=350e3J/kg),
                 Melting = MeltingParam_Caricchi()),
@@ -166,7 +170,7 @@ using WriteVTK
     MatParam                =   (
             SetMaterialParams(Name="Rock", Phase=1, 
                  Density    = ConstantDensity(ρ=2800kg/m^3),               
-               HeatCapacity = ConstantHeatCapacity(cp=1050J/kg/K),
+               HeatCapacity = ConstantHeatCapacity(Cp=1050J/kg/K),
                Conductivity = ConstantConductivity(k=1.5Watt/K/m),       
                  LatentHeat = ConstantLatentHeat(Q_L=350e3J/kg),
                     Melting = MeltingParam_Caricchi()),
diff --git a/examples/Example2D.jl b/examples/Example2D.jl
index 85dee14..f38aabc 100644
--- a/examples/Example2D.jl
+++ b/examples/Example2D.jl
@@ -18,7 +18,7 @@ using Plots
     MatParam                =   (
             SetMaterialParams(Name="Rock", Phase=1, 
                  Density    = ConstantDensity(ρ=2800kg/m^3),               
-               HeatCapacity = ConstantHeatCapacity(cp=1050J/kg/K),
+               HeatCapacity = ConstantHeatCapacity(Cp=1050J/kg/K),
                Conductivity = ConstantConductivity(k=1.5Watt/K/m),       
                  LatentHeat = ConstantLatentHeat(Q_L=350e3J/kg),
                     Melting = MeltingParam_Caricchi()),
diff --git a/examples/Example2D_ZASSy.jl b/examples/Example2D_ZASSy.jl
index 2d32c69..905074d 100644
--- a/examples/Example2D_ZASSy.jl
+++ b/examples/Example2D_ZASSy.jl
@@ -2,11 +2,11 @@
 #  It includes comparisons with 2D simulations done by the Geneva (Gregor Weber, Luca Caricchi) & UCLA (Oscar Lovera) Tracers_SimParams
 #
 # 
-const USE_GPU=true;
+const USE_GPU=false;
 using MagmaThermoKinematics
 if USE_GPU
     environment!(:gpu, Float64, 2)      # initialize parallel stencil in 2D
-    CUDA.device!(1)                     # select the GPU you use (starts @ zero)
+   # CUDA.device!(1)                     # select the GPU you use (starts @ zero)
 else
     environment!(:cpu, Float64, 2)      # initialize parallel stencil in 2D
 end
@@ -446,7 +446,7 @@ if 1==0
                             #     Conductivity = ConstantConductivity(k=3.3Watt/K/m),          # in case we use constant k
                                   Conductivity = T_Conductivity_Whittington_parameterised(),   # T-dependent k
                                  #Conductivity = T_Conductivity_Whittington(),                 # T-dependent k
-                                  HeatCapacity = ConstantHeatCapacity(cp=1000J/kg/K),
+                                  HeatCapacity = ConstantHeatCapacity(Cp=1000J/kg/K),
                                        Melting = SmoothMelting(MeltingParam_4thOrder())),      # Marxer & Ulmer melting     
                                       # Melting = MeltingParam_Caricchi()),                     # Caricchi melting
 
@@ -665,7 +665,7 @@ if 1==1
                                  Conductivity = T_Conductivity_Whittington(),                       # T-dependent k
                                   HeatCapacity = T_HeatCapacity_Whittington(),                      # T-dependent cp
                                  # Conductivity = ConstantConductivity(k=3.3Watt/K/m),               # in case we use constant k
-                                 # HeatCapacity = ConstantHeatCapacity(cp=1000J/kg/K),
+                                 # HeatCapacity = ConstantHeatCapacity(Cp=1000J/kg/K),
                                  Melting = MeltingParam_Assimilation()                              # Quadratic parameterization as in Tierney et al.
                                  #Melting = MeltingParam_Caricchi()
                                   ),       
@@ -676,7 +676,7 @@ if 1==1
                                   Conductivity = T_Conductivity_Whittington(),                       # T-dependent k
                                   HeatCapacity = T_HeatCapacity_Whittington(),                       # T-dependent cp
                                   #Conductivity = ConstantConductivity(k=3.3Watt/K/m),               # in case we use constant k
-                                  #HeatCapacity = ConstantHeatCapacity(cp=1000J/kg/K),
+                                  #HeatCapacity = ConstantHeatCapacity(Cp=1000J/kg/K),
                                         Melting = SmoothMelting(MeltingParam_Quadratic(T_s=(700+273.15)K, T_l=(1100+273.15)K))
                                         #Melting = MeltingParam_Caricchi()
                                         )       
diff --git a/examples/Example3D.jl b/examples/Example3D.jl
index 9297acc..b3e045c 100644
--- a/examples/Example3D.jl
+++ b/examples/Example3D.jl
@@ -19,7 +19,7 @@ using WriteVTK
     MatParam                =   (
             SetMaterialParams(Name="Rock", Phase=1, 
                  Density    = ConstantDensity(ρ=2800kg/m^3),               
-               HeatCapacity = ConstantHeatCapacity(cp=1050J/kg/K),
+               HeatCapacity = ConstantHeatCapacity(Cp=1050J/kg/K),
                Conductivity = ConstantConductivity(k=1.5Watt/K/m),       
                  LatentHeat = ConstantLatentHeat(Q_L=350e3J/kg),
                     Melting = MeltingParam_Caricchi()),
diff --git a/examples/Example3D_v2.jl b/examples/Example3D_v2.jl
index 3d6b1fa..4fd8f60 100644
--- a/examples/Example3D_v2.jl
+++ b/examples/Example3D_v2.jl
@@ -19,7 +19,7 @@ using WriteVTK
     MatParam                =   (
             SetMaterialParams(Name="Rock", Phase=1, 
                  Density    = ConstantDensity(ρ=2800kg/m^3),               
-               HeatCapacity = ConstantHeatCapacity(cp=1050J/kg/K),
+               HeatCapacity = ConstantHeatCapacity(Cp=1050J/kg/K),
                Conductivity = ConstantConductivity(k=1.5Watt/K/m),       
                  LatentHeat = ConstantLatentHeat(Q_L=350e3J/kg),
                     Melting = MeltingParam_Caricchi()),
diff --git a/examples/MTK_GMG_2D_example1.jl b/examples/MTK_GMG_2D_example1.jl
index 669c200..7a42059 100644
--- a/examples/MTK_GMG_2D_example1.jl
+++ b/examples/MTK_GMG_2D_example1.jl
@@ -69,7 +69,7 @@ function MTK_GMG.MTK_initialize!(Arrays::NamedTuple, Grid::GridData, Num::Numeri
      # open pvd file if requested
      if Num.Output_VTK 
         name =  joinpath(Num.SimName,Num.SimName*".pvd")
-        Num.pvd = Movie_Paraview(name=name, Initialize=true);
+        Num.pvd = movie_paraview(name=name, Initialize=true);
     end
 
     return nothing
diff --git a/examples/MTK_GMG_2D_example2.jl b/examples/MTK_GMG_2D_example2.jl
index 7d03655..59cb646 100644
--- a/examples/MTK_GMG_2D_example2.jl
+++ b/examples/MTK_GMG_2D_example2.jl
@@ -1,10 +1,4 @@
-# This is a second example that shows how to use MTK in combination with the GeophysicalModelGenerator
-# Compared to example 1, we show a few more features:
-# - How to define a custom structure with temporal values & use it in the code
-# - How to generate a model setup using GMG
-# - How to overwrite some of the default functions to customize the simulation
-# - How to create paraview files 
-
+# Unzen setup
 const USE_GPU=false;
 using MagmaThermoKinematics
 if USE_GPU
@@ -18,33 +12,34 @@ using Plots
 using Random
 using GeophysicalModelGenerator
 
-# Test setup
-println("Example 2 of the MTK - GMG integration")
-
+# Model setup
 println(" --- Generating Setup --- ")
 
 # Topography and project it. 
 
 # NOTE: The first time you do this, please set this to true, which will download the topography data from the internet and save it in a file
 if false
- using GMT, Statistics
- Topo = ImportTopo(lon = [130.0, 130.5], lat=[32.55, 32.90], file="@earth_relief_03s.grd")
- proj =  ProjectionPoint(; Lat=mean(Topo.lat.val), Lon=mean(Topo.lon.val))
- Topo_cart = Convert2CartData(Topo, proj)
- save_GMG("examples/Topo_cart", Topo_cart)
+    using GMT, Statistics
+    Topo       =   ImportTopo(lon = [130.0, 130.5], lat=[32.55, 32.90], file="@earth_relief_03s.grd")
+    proj       =   ProjectionPoint(; Lat=mean(Topo.lat.val), Lon=mean(Topo.lon.val))
+    Topo_cart  =   Convert2CartData(Topo, proj)
+    Xt,Yt,Zt   =   xyz_grid(-23:.1:23,-19:.1:19,0)
+    Topo_cart  =   ProjectCartData(CartData(Xt,Yt,Zt,(Zt=Zt,)), Topo, proj)
+
+    save_GMG("Topo_cart", Topo_cart)
 end
 Topo_cart = load_GMG("Topo_cart")       
 
 # Create 3D grid of the region
-X,Y,Z       =   XYZGrid(-23:.1:23,-19:.1:19,-20:.1:5)
+X,Y,Z       =   xyz_grid(-23:.1:23,-19:.1:19,-20:.1:5)
 Data_set3D  =   CartData(X,Y,Z,(Phases=zeros(Int64,size(X)),Temp=zeros(size(X))));       # 3D dataset
 
 # Create 2D cross-section
 Nx      =   135*6;  # resolution in x
 Nz      =   135*4;
 Data_2D =   CrossSection(Data_set3D, Start=(-20,4), End=(20,4), dims=(Nx, Nz))
-Data_2D =   AddField(Data_2D,"FlatCrossSection", FlattenCrossSection(Data_2D))
-Data_2D =   AddField(Data_2D,"Phases", Int64.(Data_2D.fields.Phases))
+Data_2D =   addfield(Data_2D,"FlatCrossSection", FlattenCrossSection(Data_2D))
+Data_2D =   addfield(Data_2D,"Phases", Int64.(Data_2D.fields.Phases))
 
 # Intersect with topography
 Below   =   BelowSurface(Data_2D, Topo_cart)
@@ -63,7 +58,6 @@ x_c, z_c, r = -10, -15, 2.5
 Volume      = 4/3*pi*r^3 # equivalent 3D volume of the anomaly [km^3]
 @views Data_2D.fields.Temp[(Data_2D.x.val .- x_c).^2 .+ (Data_2D.z.val .- z_c).^2 .< r^2] .= 800.0
 
-
 println(" --- Performing MTK models --- ")
 
 # Overwrite some of the default functions
@@ -167,24 +161,24 @@ MatParam     = (SetMaterialParams(Name="Air", Phase=0,
                                 Density      = ConstantDensity(ρ=2700kg/m^3),
                                 LatentHeat   = ConstantLatentHeat(Q_L=0.0J/kg),
                                 Conductivity = ConstantConductivity(k=3Watt/K/m),          # in case we use constant k
-                                HeatCapacity = ConstantHeatCapacity(cp=1000J/kg/K),
+                                HeatCapacity = ConstantHeatCapacity(Cp=1000J/kg/K),
                                 Melting      = SmoothMelting(MeltingParam_4thOrder())),          # Marxer & Ulmer melting     
                 SetMaterialParams(Name="Crust", Phase=1, 
                                 Density      = ConstantDensity(ρ=2700kg/m^3),
                                 LatentHeat   = ConstantLatentHeat(Q_L=3.13e5J/kg),
                                 Conductivity = T_Conductivity_Whittington_parameterised(),   # T-dependent k
-                                HeatCapacity = ConstantHeatCapacity(cp=1000J/kg/K),
+                                HeatCapacity = ConstantHeatCapacity(Cp=1000J/kg/K),
                                 Melting      = SmoothMelting(MeltingParam_4thOrder())),      # Marxer & Ulmer melting 
                 SetMaterialParams(Name="Mantle", Phase=2, 
                                 Density      = ConstantDensity(ρ=2700kg/m^3),
                                 LatentHeat   = ConstantLatentHeat(Q_L=3.13e5J/kg),
                                 Conductivity = T_Conductivity_Whittington_parameterised(),   # T-dependent k
-                                HeatCapacity = ConstantHeatCapacity(cp=1000J/kg/K)),
+                                HeatCapacity = ConstantHeatCapacity(Cp=1000J/kg/K)),
                 SetMaterialParams(Name="Dikes", Phase=3, 
                                 Density      = ConstantDensity(ρ=2700kg/m^3),
                                 LatentHeat   = ConstantLatentHeat(Q_L=3.13e5J/kg),
                                 Conductivity = T_Conductivity_Whittington_parameterised(),   # T-dependent k
-                                HeatCapacity = ConstantHeatCapacity(cp=1000J/kg/K),
+                                HeatCapacity = ConstantHeatCapacity(Cp=1000J/kg/K),
                                 Melting      = SmoothMelting(MeltingParam_4thOrder()))           # Marxer & Ulmer melting       
                 )
 
diff --git a/src/MTK_GMG.jl b/src/MTK_GMG.jl
index 6f0c26a..85d0179 100644
--- a/src/MTK_GMG.jl
+++ b/src/MTK_GMG.jl
@@ -138,7 +138,7 @@ function MTK_initialize!(Arrays::NamedTuple, Grid::GridData, Num::NumericalParam
     # Open pvd file if requested
     if Num.Output_VTK
         name =  joinpath(Num.SimName,Num.SimName*".pvd")
-        Num.pvd = Movie_Paraview(name=name, Initialize=true);
+        Num.pvd = movie_paraview(name=name, Initialize=true);
     end
 
     return nothing
@@ -201,7 +201,7 @@ function MTK_initialize!(Arrays::NamedTuple, Grid::GridData, Num::NumericalParam
     # open pvd file if requested
     if Num.Output_VTK 
         name =  joinpath(Num.SimName,Num.SimName*".pvd")
-        Num.pvd = Movie_Paraview(name=name, Initialize=true);
+        Num.pvd = movie_paraview(name=name, Initialize=true);
     end
 
     return nothing
@@ -215,7 +215,7 @@ Finalize model run
 """
 function MTK_finalize!(Arrays::NamedTuple, Grid::GridData, Num::NumericalParameters, Tracers::StructArray, Dikes::DikeParameters, CartData_input::Union{Nothing,CartData})
     if Num.Output_VTK & !isnothing(Num.pvd)
-        Movie_Paraview(pvd=Num.pvd, Finalize=true)
+        movie_paraview(pvd=Num.pvd, Finalize=true)
     end
 
     return nothing
@@ -270,9 +270,9 @@ function MTK_save_output(Grid::GridData, Arrays::NamedTuple, Tracers::StructArra
                 Data_set3D  = CartData_input
             else
                 if length(Grid.coord1D)==3
-                    X,Y,Z   =   XYZGrid(Grid.coord1D...)
+                    X,Y,Z   =   xyz_grid(Grid.coord1D...)
                 elseif length(Grid.coord1D)==2
-                    X,Y,Z   =   XYZGrid(Grid.coord1D[1], 0, Grid.coord1D[2])
+                    X,Y,Z   =   xyz_grid(Grid.coord1D[1], 0, Grid.coord1D[2])
                 end
                 Data_set3D  =   CartData(X/1e3,Y/1e3,Z/1e3, (Z=Z,))
             end
@@ -282,7 +282,7 @@ function MTK_save_output(Grid::GridData, Arrays::NamedTuple, Tracers::StructArra
             Data_set3D = add_data_CartData(Data_set3D, "MeltFraction", Array(Arrays.ϕ));
 
             # Save output to CartData
-            Num.pvd  = Write_Paraview(Data_set3D, name, pvd=Num.pvd,time=Num.time/SecYear/1e3);
+            Num.pvd  = write_paraview(Data_set3D, name, pvd=Num.pvd,time=Num.time/SecYear/1e3);
         end
     end
     return nothing
@@ -304,7 +304,7 @@ function add_data_CartData(d::CartData, name::String, data::Array)
     else
         a = data
     end
-    d = AddField(d, name, a)
+    d = addfield(d, name, a)
     return d
 end
 
diff --git a/src/MTK_GMG_2D.jl b/src/MTK_GMG_2D.jl
index c2864f5..ec6c963 100644
--- a/src/MTK_GMG_2D.jl
+++ b/src/MTK_GMG_2D.jl
@@ -53,7 +53,7 @@ There are a few functions that you can overwrite in your user code to customize
     if !isnothing(CartData_input)
        
         if !hasfield(typeof(CartData_input.fields),:FlatCrossSection)
-           error("You should add a Field :FlatCrossSection to your data structure with Data_Cross = AddField(Data_Cross,\"FlatCrossSection\", FlattenCrossSection(Data_Cross))")
+           error("You should add a Field :FlatCrossSection to your data structure with Data_Cross = addfield(Data_Cross,\"FlatCrossSection\", FlattenCrossSection(Data_Cross))")
         end
 
         Num = MTK_GMG.Setup_Model_CartData(CartData_input, Num, Mat_tup)
diff --git a/test/test_MTK_GMG_2D.jl b/test/test_MTK_GMG_2D.jl
index 3a0abf2..c90bdc6 100644
--- a/test/test_MTK_GMG_2D.jl
+++ b/test/test_MTK_GMG_2D.jl
@@ -58,7 +58,7 @@ MatParam     = (SetMaterialParams(Name="Rock & partial melt", Phase=1,
                         #     Conductivity = ConstantConductivity(k=3.3Watt/K/m),          # in case we use constant k
                             Conductivity = T_Conductivity_Whittington_parameterised(),   # T-dependent k
                             #Conductivity = T_Conductivity_Whittington(),                 # T-dependent k
-                            HeatCapacity = ConstantHeatCapacity(cp=1000J/kg/K),
+                            HeatCapacity = ConstantHeatCapacity(Cp=1000J/kg/K),
                                 Melting = SmoothMelting(MeltingParam_4thOrder())),      # Marxer & Ulmer melting     
                                 # Melting = MeltingParam_Caricchi()),                     # Caricchi melting
                 # add more parameters here, in case you have >1 phase in the model                                    
@@ -75,15 +75,15 @@ Grid, Arrays, Tracers, Dikes, time_props = MTK_GeoParams_2D(MatParam, Num, Dike_
 Topo_cart = load_GMG("../examples/Topo_cart")       # Note: Laacher seee is around [10,20]
 
 # Create 3D grid of the region
-X,Y,Z       =   XYZGrid(-23:.1:23,-19:.1:19,-20:.1:5)
+X,Y,Z       =   xyz_grid(-23:.1:23,-19:.1:19,-20:.1:5)
 Data_set3D  =   CartData(X,Y,Z,(Phases=zeros(Int64,size(X)),Temp=zeros(size(X))));       # 3D dataset
 
 # Create 2D cross-section
 Nx          =   Num.Nx;  # resolution in x
 Nz          =   Num.Nz;
 Data_2D     =   CrossSection(Data_set3D, Start=(-20,4), End=(20,4), dims=(Nx, Nz))
-Data_2D     =   AddField(Data_2D,"FlatCrossSection", FlattenCrossSection(Data_2D))
-Data_2D     =   AddField(Data_2D,"Phases", Int64.(Data_2D.fields.Phases))
+Data_2D     =   addfield(Data_2D,"FlatCrossSection", FlattenCrossSection(Data_2D))
+Data_2D     =   addfield(Data_2D,"Phases", Int64.(Data_2D.fields.Phases))
 
 # Intersect with topography
 Below = BelowSurface(Data_2D, Topo_cart)
@@ -141,7 +141,7 @@ MatParam     = (SetMaterialParams(Name="Air", Phase=0,
                                 Density    = ConstantDensity(ρ=2700kg/m^3),
                                 LatentHeat = ConstantLatentHeat(Q_L=0.0J/kg),
                                 Conductivity = ConstantConductivity(k=3Watt/K/m),          # in case we use constant k
-                                HeatCapacity = ConstantHeatCapacity(cp=1000J/kg/K),
+                                HeatCapacity = ConstantHeatCapacity(Cp=1000J/kg/K),
                                 Melting = SmoothMelting(MeltingParam_4thOrder())),          # Marxer & Ulmer melting     
 
                 SetMaterialParams(Name="Crust", Phase=1, 
@@ -149,14 +149,14 @@ MatParam     = (SetMaterialParams(Name="Air", Phase=0,
                                 LatentHeat = ConstantLatentHeat(Q_L=3.13e5J/kg),
                                 Conductivity = T_Conductivity_Whittington_parameterised(),   # T-dependent k
                                 #Conductivity = T_Conductivity_Whittington(),                 # T-dependent k
-                                HeatCapacity = ConstantHeatCapacity(cp=1000J/kg/K),
+                                HeatCapacity = ConstantHeatCapacity(Cp=1000J/kg/K),
                                 Melting = SmoothMelting(MeltingParam_4thOrder())),      # Marxer & Ulmer melting 
 
                 SetMaterialParams(Name="Mantle", Phase=2, 
                                 Density    = ConstantDensity(ρ=2700kg/m^3),
                                 LatentHeat = ConstantLatentHeat(Q_L=3.13e5J/kg),
                                 Conductivity = T_Conductivity_Whittington_parameterised(),   # T-dependent k
-                                HeatCapacity = ConstantHeatCapacity(cp=1000J/kg/K)),
+                                HeatCapacity = ConstantHeatCapacity(Cp=1000J/kg/K)),
 
                 SetMaterialParams(Name="Dikes", Phase=3, 
                                 Density    = ConstantDensity(ρ=2700kg/m^3),
@@ -164,7 +164,7 @@ MatParam     = (SetMaterialParams(Name="Air", Phase=0,
                         #     Conductivity = ConstantConductivity(k=3.3Watt/K/m),          # in case we use constant k
                                 Conductivity = T_Conductivity_Whittington_parameterised(),   # T-dependent k
                                 #Conductivity = T_Conductivity_Whittington(),                 # T-dependent k
-                                HeatCapacity = ConstantHeatCapacity(cp=1000J/kg/K),
+                                HeatCapacity = ConstantHeatCapacity(Cp=1000J/kg/K),
                                 Melting = SmoothMelting(MeltingParam_4thOrder()))      # Marxer & Ulmer melting     
                                   
                 )
diff --git a/test/test_MTK_GMG_3D.jl b/test/test_MTK_GMG_3D.jl
index c66d5a4..b17b849 100644
--- a/test/test_MTK_GMG_3D.jl
+++ b/test/test_MTK_GMG_3D.jl
@@ -66,7 +66,7 @@ MatParam     = (SetMaterialParams(Name="Rock & partial melt", Phase=1,
                         #     Conductivity = ConstantConductivity(k=3.3Watt/K/m),          # in case we use constant k
                             Conductivity = T_Conductivity_Whittington_parameterised(),   # T-dependent k
                             #Conductivity = T_Conductivity_Whittington(),                 # T-dependent k
-                            HeatCapacity = ConstantHeatCapacity(cp=1000J/kg/K),
+                            HeatCapacity = ConstantHeatCapacity(Cp=1000J/kg/K),
                                 Melting = SmoothMelting(MeltingParam_4thOrder())),      # Marxer & Ulmer melting     
                                 # Melting = MeltingParam_Caricchi()),                     # Caricchi melting
                 # add more parameters here, in case you have >1 phase in the model                                    
@@ -84,7 +84,7 @@ Topo_cart = load_GMG("../examples/Topo_cart")       # Note: Laacher seee is arou
 
 # Create 3D grid of the region
 Nx,Ny,Nz = 100,100,100
-X,Y,Z       =   XYZGrid(range(-23,23, length=Nx),range(-19,19, length=Ny),range(-20,5, length=Nz))
+X,Y,Z       =   xyz_grid(range(-23,23, length=Nx),range(-19,19, length=Ny),range(-20,5, length=Nz))
 Data_3D     =   CartData(X,Y,Z,(Phases=zeros(Int64,size(X)),Temp=zeros(size(X))));       # 3D dataset
 
 # Intersect with topography
@@ -131,7 +131,7 @@ MatParam     = (SetMaterialParams(Name="Air", Phase=0,
                                 Density    = ConstantDensity(ρ=2700kg/m^3),
                                 LatentHeat = ConstantLatentHeat(Q_L=0.0J/kg),
                                 Conductivity = ConstantConductivity(k=3Watt/K/m),          # in case we use constant k
-                                HeatCapacity = ConstantHeatCapacity(cp=1000J/kg/K),
+                                HeatCapacity = ConstantHeatCapacity(Cp=1000J/kg/K),
                                 Melting = SmoothMelting(MeltingParam_4thOrder())),          # Marxer & Ulmer melting     
 
                 SetMaterialParams(Name="Crust", Phase=1, 
@@ -139,14 +139,14 @@ MatParam     = (SetMaterialParams(Name="Air", Phase=0,
                                 LatentHeat = ConstantLatentHeat(Q_L=3.13e5J/kg),
                                 Conductivity = T_Conductivity_Whittington_parameterised(),   # T-dependent k
                                 #Conductivity = T_Conductivity_Whittington(),                 # T-dependent k
-                                HeatCapacity = ConstantHeatCapacity(cp=1000J/kg/K),
+                                HeatCapacity = ConstantHeatCapacity(Cp=1000J/kg/K),
                                 Melting = SmoothMelting(MeltingParam_4thOrder())),      # Marxer & Ulmer melting 
 
                 SetMaterialParams(Name="Mantle", Phase=2, 
                                 Density    = ConstantDensity(ρ=2700kg/m^3),
                                 LatentHeat = ConstantLatentHeat(Q_L=3.13e5J/kg),
                                 Conductivity = T_Conductivity_Whittington_parameterised(),   # T-dependent k
-                                HeatCapacity = ConstantHeatCapacity(cp=1000J/kg/K)),
+                                HeatCapacity = ConstantHeatCapacity(Cp=1000J/kg/K)),
 
                 SetMaterialParams(Name="Dikes", Phase=3, 
                                 Density    = ConstantDensity(ρ=2700kg/m^3),
@@ -154,7 +154,7 @@ MatParam     = (SetMaterialParams(Name="Air", Phase=0,
                         #     Conductivity = ConstantConductivity(k=3.3Watt/K/m),          # in case we use constant k
                                 Conductivity = T_Conductivity_Whittington_parameterised(),   # T-dependent k
                                 #Conductivity = T_Conductivity_Whittington(),                 # T-dependent k
-                                HeatCapacity = ConstantHeatCapacity(cp=1000J/kg/K),
+                                HeatCapacity = ConstantHeatCapacity(Cp=1000J/kg/K),
                                 Melting = SmoothMelting(MeltingParam_4thOrder()))      # Marxer & Ulmer melting     
                                   
                 )
diff --git a/test/test_ZASSY.jl b/test/test_ZASSY.jl
index 9501255..616e558 100644
--- a/test/test_ZASSY.jl
+++ b/test/test_ZASSY.jl
@@ -355,7 +355,7 @@ if 1==1
                             #     Conductivity = ConstantConductivity(k=3.3Watt/K/m),          # in case we use constant k
                                 Conductivity = T_Conductivity_Whittington_parameterised(),   # T-dependent k
                                 #Conductivity = T_Conductivity_Whittington(),                 # T-dependent k
-                                HeatCapacity = ConstantHeatCapacity(cp=1000J/kg/K),
+                                HeatCapacity = ConstantHeatCapacity(Cp=1000J/kg/K),
                                     Melting = SmoothMelting(MeltingParam_4thOrder())),      # Marxer & Ulmer melting     
                                     # Melting = MeltingParam_Caricchi()),                     # Caricchi melting
                     # add more parameters here, in case you have >1 phase in the model