diff --git a/Project.toml b/Project.toml
index 237aabbb..a306dc47 100644
--- a/Project.toml
+++ b/Project.toml
@@ -16,11 +16,13 @@ Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
 IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
+NetCDF = "30363a11-5582-574a-97bb-aa9a979735b9"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 OnlineStats = "a15396b6-48d5-5d58-9928-6d29437db91e"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 ParallelUtilities = "fad6cfc8-4f83-11e9-06cc-151124046ad0"
+PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
@@ -29,6 +31,7 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 WeightedOnlineStats = "bbac0a1f-7c9d-5672-960b-c6ca726e5d5d"
 YAXArrayBase = "90b8fcef-0c2d-428d-9c56-5f86629e9d14"
+Zarr = "0a941bbe-ad1d-11e8-39d9-ab76183a1d99"
 
 [compat]
 CFTime = "0.0, 0.1"
diff --git a/src/Precompile/precompile.jl b/src/Precompile/precompile.jl
new file mode 100644
index 00000000..f3633f45
--- /dev/null
+++ b/src/Precompile/precompile.jl
@@ -0,0 +1,55 @@
+using PrecompileTools
+
+@recompile_invalidations begin
+    # load packages
+    using Dates
+    using Statistics
+    using Zarr
+end
+
+
+@setup_workload begin
+    # Putting some things in `@setup_workload` instead of `@compile_workload` can reduce the size of the
+    # precompile file and potentially make loading faster.      
+    @compile_workload begin
+        z0 = YAXArray(rand(10,10))
+        nameout = string(tempname(),".zarr")
+        z0 = YAXArray(rand(10,10))
+        z1 = zcreate(Float64, 10,10, path=nameout)
+        z1 .= z0[:,:]
+        
+        
+        
+        
+        #savecube(z0, nameout)
+        # c = Cube(nameout)
+        # show(c[1,1])  
+        # show(c[1,1:3])
+        # show(c[2,1:3])
+        # show(c)
+        # csum = mapslices(sum, c, dims="Dim_1")
+        # cmean = mapslices(mean, c, dims="Dim_1")
+       
+       
+        # z0 = YAXArray(rand(10,10))
+        # nameout = string(tempname(),".zarr")
+        # z1 = zcreate(Float64, 10,10, path=nameout)
+        # z1 .= z0[:,:]
+        # z2 = zopen(nameout)
+        # show(z2[1,1])  
+        # show(z2[1,1:3])
+        # show(z2[2,1:3])
+        # show(z2)
+        # zsum = mapslices(sum, z2, dims=1)
+        # zmean = mapslices(mean, z2, dims=2)
+
+
+        # a = YAXArray(rand(10,10))
+        # show(a[1,1])  
+        # show(a[Dim_1=1..3])
+        # show(a[Dim_2=1..3])
+        # show(a)
+        # a1 = mapslices(sum, a, dims="Dim_1")
+        # a3 = mapslices(mean, a, dims="Dim_2")
+    end
+end
diff --git a/src/Precompile/precompile_zarr.jl b/src/Precompile/precompile_zarr.jl
new file mode 100644
index 00000000..40d6a4ce
--- /dev/null
+++ b/src/Precompile/precompile_zarr.jl
@@ -0,0 +1,45 @@
+using Pkg
+
+using PrecompileTools
+
+@recompile_invalidations begin
+    # load packages
+    using Dates
+    using Statistics
+    using Zarr
+    # using YAXArrays
+end
+
+
+@setup_workload begin
+    # Putting some things in `@setup_workload` instead of `@compile_workload` can reduce the size of the
+    # precompile file and potentially make loading faster.      
+    @compile_workload begin
+
+        t =  Date("2020-01-01"):Month(1):Date("2022-12-31")
+
+        ## create cube axes
+        axes = [RangeAxis("Lon", 1:10), RangeAxis("Lat", 1:10), RangeAxis("Time", t)]
+
+        ## assign values to a cube
+        z0 = YAXArray(axes, reshape(1:3600, (10,10,36)))
+
+        # saving and opening Zarr files
+        nameout = string(tempname(),".zarr")
+        z1 = zcreate(Float64, 10,10,36, path=nameout)
+        z1 .= z0[:,:,:]
+        z2 = zopen(nameout)
+
+        # simple data visualization
+        show(z2[1,1,1])
+        show(z2[1,1:3,1])
+        show(z2[2,1:3,1])
+        show(z2)
+
+        # simple statistics across different axes
+        zsum = mapslices(sum, z2, dims=1)
+        zmean = mapslices(mean, z2, dims=2)
+        zmean = mapslices(mean, z2, dims=3)
+
+    end 
+end
\ No newline at end of file
diff --git a/src/YAXArrays.jl b/src/YAXArrays.jl
index cf24d616..4af096b6 100644
--- a/src/YAXArrays.jl
+++ b/src/YAXArrays.jl
@@ -1,4 +1,5 @@
 module YAXArrays
+using Zarr
 
 """
 Default configuration for YAXArrays, has the following fields:
@@ -50,7 +51,7 @@ using YAXArrayBase: getattributes
 @reexport using .Datasets
 # from YAXTools
 
-# include("precompile.jl")
+include("Precompile/precompile.jl")
 # _precompile_()
 
 end # module