diff --git a/Project.toml b/Project.toml
index 2afebacf..85d31495 100644
--- a/Project.toml
+++ b/Project.toml
@@ -52,7 +52,7 @@ FFMPEG = "0.4"
 FileIO = "1"
 GDAL_jll = "300 - 304"
 GLMakie = "0.8 - 0.11, 0.13"
-GMT = "1"
+GMT = ">=1.17"
 GeoParams = "0.2 - 0.7"
 Geodesy = "1"
 GeometryBasics = "0.1 - 0.5"
diff --git a/src/LaMEM_io.jl b/src/LaMEM_io.jl
index d3ba8577..c6ee316a 100644
--- a/src/LaMEM_io.jl
+++ b/src/LaMEM_io.jl
@@ -406,19 +406,28 @@ function save_LaMEM_markers_parallel(Grid::CartData; PartitioningFile = empty, d
     y = ustrip.(Grid.y.val[1, :, 1])
     z = ustrip.(Grid.z.val[1, 1, :])
 
-    if haskey(Grid.fields, :Phases)
-        Phases = Grid.fields[:Phases]
+    Phases = if haskey(Grid.fields, :Phases)
+        Grid.fields[:Phases]
     else
         error("You must provide the field :Phases in the structure")
     end
 
-    if haskey(Grid.fields, :Temp)
-        Temp = Grid.fields[:Temp]
+    Temp = if haskey(Grid.fields, :Temp)
+        Grid.fields[:Temp]
     else
         if verbose
             println("Field :Temp is not provided; setting it to zero")
         end
-        Temp = zeros(size(Phases))
+        zeros(size(Phases))
+    end
+
+    APS = if haskey(Grid.fields, :APS)
+        Grid.fields[:APS]
+    else
+        if verbose
+            println("Field :APS is not provided; setting it to zero")
+        end
+        zeros(size(Phases))
     end
 
     if PartitioningFile == empty
@@ -459,10 +468,11 @@ function save_LaMEM_markers_parallel(Grid::CartData; PartitioningFile = empty, d
         part_z = ustrip.(Grid.z.val[x_start[n]:x_end[n], y_start[n]:y_end[n], z_start[n]:z_end[n]])
         part_phs = Phases[x_start[n]:x_end[n], y_start[n]:y_end[n], z_start[n]:z_end[n]]
         part_T = Temp[x_start[n]:x_end[n], y_start[n]:y_end[n], z_start[n]:z_end[n]]
+        part_APS = APS[x_start[n]:x_end[n], y_start[n]:y_end[n], z_start[n]:z_end[n]]
         num_particles = size(part_x, 1) * size(part_x, 2) * size(part_x, 3)
 
         # Information vector per processor
-        num_prop = 5       # number of properties we save [x/y/z/phase/T]
+        num_prop = 6       # number of properties we save [x/y/z/phase/T/APS]
         lvec_info = num_particles
 
         lvec_prtcls = zeros(Float64, num_prop * num_particles)
@@ -472,6 +482,7 @@ function save_LaMEM_markers_parallel(Grid::CartData; PartitioningFile = empty, d
         lvec_prtcls[3:num_prop:end] = part_z[:]
         lvec_prtcls[4:num_prop:end] = part_phs[:]
         lvec_prtcls[5:num_prop:end] = part_T[:]
+        lvec_prtcls[6:num_prop:end] = part_APS[:]
 
         # Write output files
         if ~isdir(directory)
@@ -552,7 +563,10 @@ function PetscBinaryWrite_Vec(filename, A)
         n = length(A)
         nummark = A[1]            # number of markers
 
-        write(f, hton(Float64(1211214)))     # header (not actually used)
+        # header number encodes the version of particle data file
+        # version with APS is 1211215
+        # version without APS was 1211214
+        write(f, hton(Float64(1211215)))     # header
         write(f, hton(Float64(nummark)))     # info about # of markers written
 
         for i in 2:n