diff --git a/Project.toml b/Project.toml
index 1576262d..ad4c4201 100644
--- a/Project.toml
+++ b/Project.toml
@@ -17,12 +17,14 @@ OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [weakdeps]
-Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 Climatology = "9e9a4d37-2d2e-41e3-8b85-f7978328d9c7"
+Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+MITgcm = "dce5fa8e-68ce-4431-a242-9469c69627a0"
 
 [extensions]
-IndividualDisplacementsMakieExt = ["Makie"]
 IndividualDisplacementsClimatologyExt = ["Climatology"]
+IndividualDisplacementsMakieExt = ["Makie"]
+IndividualDisplacementsMITgcmExt = ["MITgcm"]
 
 [compat]
 CSV = "0.10"
@@ -35,6 +37,7 @@ Glob = "1"
 JLD2 = "0.5"
 Makie = "0.21"
 MeshArrays = "0.3"
+MITgcm = "0.5"
 NetCDF = "0.12"
 OrdinaryDiffEq = "6"
 julia = "^1.10"
diff --git a/examples/worldwide/global_ocean_circulation.jl b/examples/worldwide/global_ocean_circulation.jl
index d6cf1bf7..86692933 100644
--- a/examples/worldwide/global_ocean_circulation.jl
+++ b/examples/worldwide/global_ocean_circulation.jl
@@ -19,12 +19,21 @@ end
 # ╔═╡ 1821a034-5a3f-4a3a-a971-a80a26aa01cb
 using Pkg; Pkg.status()
 
-# ╔═╡ 104ce9b0-3fd1-11ec-3eff-3b029552e3d9
+# ╔═╡ 2acbfc79-3926-4c5a-9994-e3222c60c377
 begin
-	using Statistics, PlutoUI
+	import IndividualDisplacements, Climatology, MITgcm
+    using Statistics, PlutoUI, CairoMakie
+
+    ECCOmodule = IndividualDisplacements.ECCO
+    CSV = IndividualDisplacements.CSV
+    DataFrames = IndividualDisplacements.DataFrames
+
 	output_path=joinpath(tempdir(),"global_ocean_tmp")
 	!isdir(output_path) ? mkdir(output_path) : nothing
-end
+
+    path_input = IndividualDisplacements.datadeps.getdata("global_ocean_circulation_inputs")
+    path_replay = IndividualDisplacements.datadeps.getdata("global_ocean_circulation_outputs")
+end	
 
 # ╔═╡ c9e9faa8-f5f0-479c-bc85-877ff7114883
 md"""# Global Climatology
@@ -43,94 +52,6 @@ For additional documentation see :
 # ╔═╡ 171fa252-7a35-4d4a-a940-60de77327cf4
 TableOfContents()
 
-# ╔═╡ 94ca10ae-6a8a-4038-ace0-07d7d9026712
-md"""## 2. Configure `FlowFields` 
-
-####
-
-`global_ocean_circulation` reads in the global Ocean model grid and flow fields. It returns the `FlowFields` data structure `𝑃`. Ancillary variables are stored in NamedTuple `𝐷`.
-
-!!! note
-    The global grid is handled via `MeshArrays.jl` and monthly climatology fields accessed via `Climatology.jl`.
-
-```
-zrng=ARGS[1]
-zrng=parse.(Int,split(zrng,":"))
-zrng=zrng[1]:zrng[2]
-println(zrng)
-```
-
-!!! note
-    For non-interactive, one can add here code as shown above. Then `julia --project=. main.jl 21:27 &`.
-"""
-
-# ╔═╡ 2fb4d76d-56d5-4c3e-bc7c-ba0cac57e0e3
-output_path
-
-# ╔═╡ f1215951-2eb2-490b-875a-91c1205b8f63
-md"""## 3. Trajectory Computation
-
-### 3.1 Initialization
-
-- initialize individual positions (`init_gulf_stream`)
-- initial integration from time 0 to 0.5 month
-"""
-
-# ╔═╡ 6158a5e4-89e0-4496-ab4a-044d1e3e8cc0
-md""" ### 3.2 Monthly Step Function
-
-Since the Ocean circulation varies from month to month, the `𝐼.𝐷.🔄` function is used to update flow fields from one month to the next. 
-
-Time variable flow fields are easily handled by defining a `step!` function that wraps around `∫!`. Here it does three things
-
-- `𝐼.𝐷.🔄` resets the velocity inputs (`𝐼.𝐷.u0` etc) to bracket time `t_ϵ=𝐼.𝑃.𝑇[2]+eps(𝐼.𝑃.𝑇[2])`
-- `reset_📌!` randomly selects a fraction of the individuals and resets their positions. This can be useful to maintain homogeneous coverage of the domain by the fleet of individuals.
-- `∫!(𝐼)` solves for the individual trajectories over one month (`𝐼.𝑃.𝑇`) with updated velocity fields (`𝐼.𝑃.u0` etc). It also adds diagnostics to the `DataFrame` used to record variables along the trajectory (`𝐼.🔴`).
-"""
-
-# ╔═╡ 7efadea7-4542-40cf-893a-40a75e9c52be
-md"""### 3.3 Monthly Simulation Loop
-
-!!! note
-    `add_lonlat!` derives geographic locations (longitude and latitude) from local grid coordinates (x, y, etc).
-"""
-
-# ╔═╡ c5ba37e9-2a68-4448-a2cb-dea1fbf08f1e
-md"""## 4. Visualize Displacements
-
-!!! note
-    For offline visualization, see code below.
-"""
-
-# ╔═╡ 1a6af0eb-ab2a-4999-8063-f218b2f3f651
-"output path is $(output_path)"
-
-# ╔═╡ 15077957-64d5-46a5-8a87-a76ad619cf38
-md"""## 5. Summary Statistics
-
-Here we briefly demontrate the use of [DataFrames.jl](https://juliadata.github.io/DataFrames.jl/latest/) to analyze the output (𝐼.🔴) of our simulation.
-"""
-
-# ╔═╡ 6bcd1a14-6bee-4549-8007-61952909b18f
-md"""## Appendix"""
-
-# ╔═╡ 2acbfc79-3926-4c5a-9994-e3222c60c377
-module inc 
-	import IndividualDisplacements
-	import Climatology, MITgcm, MeshArrays
-    import CairoMakie, Makie
-    import DataFrames, CSV, JLD2
-	import FileIO, NetCDF, DataDeps, Colors
-
-    p0=joinpath(dirname(pathof(IndividualDisplacements)),"..","examples")
-    f0=joinpath(p0,"worldwide","ECCO_FlowFields.jl")
-    f1=( isfile("ECCO_FlowFields.jl") ? "ECCO_FlowFields.jl" :  f0 )
-    include(f1)
-
-    path_input = inc.IndividualDisplacements.datadeps.getdata("global_ocean_circulation_inputs")
-    path_replay = inc.IndividualDisplacements.datadeps.getdata("global_ocean_circulation_outputs")
-end	
-
 # ╔═╡ 7fec71b4-849f-4369-bec2-26bfe2e00a97
 begin
 	✔0="selected parameters"
@@ -143,7 +64,7 @@ begin
 	bind_zoomin = (@bind zoom_in CheckBox(default=false))
 	bind_zrng = @bind zrng Select([0:11,11:21,21:27],default=0:11)
 	
-    path_IC = inc.path_input
+    path_IC = path_input
     file_IC = joinpath(path_IC,"initial_10_1.csv")
 	backward_time = false
 	file_base = basename(file_IC)[1:end-4]
@@ -165,14 +86,26 @@ begin
 	"""
 end
 
-# ╔═╡ 63b68e72-76c1-4104-bf76-dd9eefc4e225
-md"""### 3.4 Replay previous simulation
+# ╔═╡ 94ca10ae-6a8a-4038-ace0-07d7d9026712
+md"""## 2. Configure `FlowFields` 
 
-If the replay option ($(bind_replay)) has been selected then we reload the result of a previous computation from file.
-"""
+####
 
-# ╔═╡ e49948eb-8b62-4579-a80c-b07624da6f3b
-md"""latitude box : $(bind_j)"""
+`global_ocean_circulation` reads in the global Ocean model grid and flow fields. It returns the `FlowFields` data structure `𝑃`. Ancillary variables are stored in NamedTuple `𝐷`.
+
+!!! note
+    The global grid is handled via `MeshArrays.jl` and monthly climatology fields accessed via `Climatology.jl`.
+
+```
+zrng=ARGS[1]
+zrng=parse.(Int,split(zrng,":"))
+zrng=zrng[1]:zrng[2]
+println(zrng)
+```
+
+!!! note
+    For non-interactive, one can add here code as shown above. Then `julia --project=. main.jl 21:27 &`.
+"""
 
 # ╔═╡ 218b9beb-68f2-4498-a96d-08e0719b4cff
 begin
@@ -182,9 +115,9 @@ begin
 	k=parse(Int,ktxt) #vertical level or 0
 	k!==0 ? zs=[k-0.5:k-0.5] : zs=zrng
 	
-	inc.Climatology.get_ecco_velocity_if_needed()
+	Climatology.get_ecco_velocity_if_needed()
 	
-	𝑃,𝐷=inc.ECCO_FlowFields.init_FlowFields(k=k,backward_time=backward_time)
+	𝑃,𝐷=ECCOmodule.init_FlowFields(k=k,backward_time=backward_time)
 
 	println("Done with Setting Up FlowFields")
 	tmp1="  - flow field depth level = "*(k==0 ? "three-dimensional" : "level $(k)") 
@@ -193,6 +126,18 @@ begin
 	println(tmp1)
 end
 
+# ╔═╡ 2fb4d76d-56d5-4c3e-bc7c-ba0cac57e0e3
+output_path
+
+# ╔═╡ f1215951-2eb2-490b-875a-91c1205b8f63
+md"""## 3. Trajectory Computation
+
+### 3.1 Initialization
+
+- initialize individual positions (`init_gulf_stream`)
+- initial integration from time 0 to 0.5 month
+"""
+
 # ╔═╡ f727992f-b72a-45bc-93f1-cc8daf89af0f
 begin
 	✔0
@@ -206,23 +151,51 @@ begin
 		nm=12 #number of months
 	end
 	
-	df = inc.IndividualDisplacements.init.init_gulf_stream(np , 𝐷, zs=zs)
+	df = IndividualDisplacements.init.init_gulf_stream(np , 𝐷, zs=zs)
 
 	if !(k==0)
-		𝑆 = inc.ECCO_FlowFields.init_storage(np,100,1,50)
-		𝐼 = Individuals(𝑃,df.x,df.y,df.fid,(𝐷=merge(𝐷,𝑆),∫=inc.ECCO_FlowFields.custom∫))
+		𝑆 = ECCOmodule.init_storage(np,100,1,50)
+		𝐼 = Individuals(𝑃,df.x,df.y,df.fid,(𝐷=merge(𝐷,𝑆),∫=ECCOmodule.custom∫))
 		my∫! = ∫!
 	else		
-		𝑆 = inc.ECCO_FlowFields.init_storage(np,100,length(𝐷.Γ.RC),50)
-		𝐼 = inc.IndividualDisplacements.Individuals(𝑃,df.x,df.y,df.z,df.fid,
-			(𝐷=merge(𝐷,𝑆),∫=inc.ECCO_FlowFields.custom∫,🔧=inc.ECCO_FlowFields.custom🔧,🔴=deepcopy(inc.ECCO_FlowFields.custom🔴)))
-		my∫! = inc.ECCO_FlowFields.custom∫!
+		𝑆 = ECCOmodule.init_storage(np,100,length(𝐷.Γ.RC),50)
+		𝐼 = Individuals(𝑃,df.x,df.y,df.z,df.fid,
+			(𝐷=merge(𝐷,𝑆),∫=ECCOmodule.custom∫,
+            🔧=ECCOmodule.custom🔧,
+            🔴=deepcopy(ECCOmodule.custom🔴)))
+		my∫! = ECCOmodule.custom∫!
 	end
 
 	📌_reference=deepcopy(𝐼.📌)
 	𝐼
 end
 
+# ╔═╡ 6158a5e4-89e0-4496-ab4a-044d1e3e8cc0
+md""" ### 3.2 Monthly Step Function
+
+Since the Ocean circulation varies from month to month, the `𝐼.𝐷.🔄` function is used to update flow fields from one month to the next. 
+
+Time variable flow fields are easily handled by defining a `step!` function that wraps around `∫!`. Here it does three things
+
+- `𝐼.𝐷.🔄` resets the velocity inputs (`𝐼.𝐷.u0` etc) to bracket time `t_ϵ=𝐼.𝑃.𝑇[2]+eps(𝐼.𝑃.𝑇[2])`
+- `reset_📌!` randomly selects a fraction of the individuals and resets their positions. This can be useful to maintain homogeneous coverage of the domain by the fleet of individuals.
+- `∫!(𝐼)` solves for the individual trajectories over one month (`𝐼.𝑃.𝑇`) with updated velocity fields (`𝐼.𝑃.u0` etc). It also adds diagnostics to the `DataFrame` used to record variables along the trajectory (`𝐼.🔴`).
+"""
+
+# ╔═╡ a2375720-f599-43b9-a7fb-af17956309b6
+function step!(𝐼::Individuals)
+    𝐼.𝐷.🔄(𝐼)
+    𝐼.𝐷.frac > 0 ? ECCOmodule.reset_📌!(𝐼,𝐼.𝐷.frac,📌_reference) : nothing
+    my∫!(𝐼)
+end
+
+# ╔═╡ 7efadea7-4542-40cf-893a-40a75e9c52be
+md"""### 3.3 Monthly Simulation Loop
+
+!!! note
+    `add_lonlat!` derives geographic locations (longitude and latitude) from local grid coordinates (x, y, etc).
+"""
+
 # ╔═╡ a3e45927-5d53-42be-b7b7-489d6e7a6fe5
 if !do_replay
 	𝑇=(0.0,𝐼.𝑃.𝑇[2])
@@ -232,13 +205,6 @@ else
 	✔1="Skipping Initial Integration (replay instead)"
 end
 
-# ╔═╡ a2375720-f599-43b9-a7fb-af17956309b6
-function step!(𝐼::inc.IndividualDisplacements.Individuals)
-    𝐼.𝐷.🔄(𝐼)
-    𝐼.𝐷.frac > 0 ? inc.ECCO_FlowFields.reset_📌!(𝐼,𝐼.𝐷.frac,📌_reference) : nothing
-    my∫!(𝐼)
-end
-
 # ╔═╡ 1044c5aa-1a56-45b6-a4c6-63d24eea878d
 if !do_replay
 	✔1	
@@ -254,7 +220,7 @@ if !do_replay
 	✔2
 	!isdir(output_path) ? mkdir(output_path) : nothing
 	file_output=joinpath(output_path,file_base*".csv")
-	inc.CSV.write(file_output, Float32.(𝐼.🔴))
+	CSV.write(file_output, Float32.(𝐼.🔴))
 	✔3="Done with saving trajectories"
 else
 	✔3="Skipping save (replay instead)"
@@ -264,25 +230,41 @@ end
 if !do_replay
 	✔3
 	file_output_csv=joinpath(output_path,file_base*".csv")
-	inc.CSV.write(file_output_csv, Float32.(𝐼.🔴))
+	CSV.write(file_output_csv, Float32.(𝐼.🔴))
 else
 	✔3
 	"Skipping File Output (replay instead)"
 end
 
+# ╔═╡ 63b68e72-76c1-4104-bf76-dd9eefc4e225
+md"""### 3.4 Replay previous simulation
+
+If the replay option ($(bind_replay)) has been selected then we reload the result of a previous computation from file.
+"""
+
 # ╔═╡ 397e5491-56ce-44ba-81d4-2982b0c3f503
 begin
 	#@bind fil_replay FilePicker()    
-    fil_replay=joinpath(inc.path_replay,"initial_10_$(j)_▶▶.csv")
+    fil_replay=joinpath(path_replay,"initial_10_$(j)_▶▶.csv")
 end
 
+# ╔═╡ c5ba37e9-2a68-4448-a2cb-dea1fbf08f1e
+md"""## 4. Visualize Displacements
+
+!!! note
+    For offline visualization, see code below.
+"""
+
+# ╔═╡ 1a6af0eb-ab2a-4999-8063-f218b2f3f651
+"output path is $(output_path)"
+
 # ╔═╡ 33fb4a15-b5ef-46f8-9e4e-c20da4536195
 if do_replay&&isa(fil_replay,Dict)&&haskey(fil_replay,"name")
 	#tmp_🔴=CSV.read(fil_replay["name"],DataFrame)
-	tmp_🔴=UInt8.(fil_replay["data"]) |> IOBuffer |> inc.CSV.File |> inc.DataFrames.DataFrame
+	tmp_🔴=UInt8.(fil_replay["data"]) |> IOBuffer |> CSV.File |> DataFrames.DataFrame
 	tmp_file_base=split(fil_replay["name"],'.')[1]
 elseif do_replay&&isa(fil_replay,String)
-	tmp_🔴=inc.CSV.read(fil_replay,inc.DataFrames.DataFrame)
+	tmp_🔴=CSV.read(fil_replay,DataFrames.DataFrame)
 #	tmp_file_base=basename(fil_replay)[1:end-4]
 	tmp_file_base=basename(fil_replay)
 	tmp_file_base=split(tmp_file_base,"_▶▶")[1]
@@ -304,14 +286,14 @@ if !isempty(tmp_🔴)
 		ylims=(20.0,67.0)
 	end
 
-    x=inc.IndividualDisplacements.InDiPlot( data=(I=𝐼,df=tmp_🔴,), options=(plot_type=:global_plot1,) )
-	fig,tt=inc.Makie.plot(x)
+    x=IndividualDisplacements.InDiPlot( data=(I=𝐼,df=tmp_🔴,), options=(plot_type=:global_plot1,) )
+	fig,tt=Makie.plot(x)
 
 	file_output_png=joinpath(output_path,tmp_file_base*".png")
-	inc.Makie.save(file_output_png,fig)
+	Makie.save(file_output_png,fig)
 
 #	file_output_mp4=joinpath(output_path,tmp_file_base*".mp4")
-#	inc.record(fig, file_output_mp4, 1:nt, framerate = 30) do t
+#	record(fig, file_output_mp4, 1:nt, framerate = 30) do t
 #		tt[]=t
 #	end
 	
@@ -320,14 +302,30 @@ else
 	"nothing to plot"
 end
 
-# ╔═╡ 0251b905-82e1-41c7-9917-dfdb980eef4f
-tmp_🔴
+# ╔═╡ e49948eb-8b62-4579-a80c-b07624da6f3b
+md"""latitude box : $(bind_j)"""
+
+# ╔═╡ 15077957-64d5-46a5-8a87-a76ad619cf38
+md"""## 5. Summary Statistics
+
+Here we briefly demontrate the use of [DataFrames.jl](https://juliadata.github.io/DataFrames.jl/latest/) to analyze the output (𝐼.🔴) of our simulation.
+"""
 
 # ╔═╡ 6e43a2af-bf01-4f42-a4ba-1874a8cf4885
 begin
 	✔2
-	gdf = inc.DataFrames.groupby(tmp_🔴, :ID)
-	sgdf= inc.DataFrames.combine(gdf,inc.DataFrames.nrow,:lat => mean)
+	gdf = DataFrames.groupby(tmp_🔴, :ID)
+	sgdf= DataFrames.combine(gdf,DataFrames.nrow,:lat => mean)
+end
+
+# ╔═╡ 0251b905-82e1-41c7-9917-dfdb980eef4f
+tmp_🔴
+
+# ╔═╡ 6bcd1a14-6bee-4549-8007-61952909b18f
+md"""## Appendix"""
+
+# ╔═╡ 104ce9b0-3fd1-11ec-3eff-3b029552e3d9
+begin
 end
 
 # ╔═╡ e9862707-3c00-434b-9501-e590e331b0b3
diff --git a/examples/worldwide/three_dimensional_ocean.jl b/examples/worldwide/three_dimensional_ocean.jl
index f73f6cd6..0b707473 100644
--- a/examples/worldwide/three_dimensional_ocean.jl
+++ b/examples/worldwide/three_dimensional_ocean.jl
@@ -6,13 +6,8 @@ using InteractiveUtils
 
 # ╔═╡ 16eab80b-325b-43bd-8bda-6b9ed27513a8
 begin
-	using IndividualDisplacements, CairoMakie, Climatology, NetCDF
-
-    p0=joinpath(dirname(pathof(IndividualDisplacements)),"..","examples")
-    f0=joinpath(p0,"worldwide","OCCA_FlowFields.jl")
-    f1=( isfile("OCCA_FlowFields.jl") ? "OCCA_FlowFields.jl" :  f0 )
-    include(f1)
-
+	using IndividualDisplacements, CairoMakie, Climatology
+    OCCAmodule=IndividualDisplacements.OCCA
 	"Done with loading packages"
 end
 
@@ -29,7 +24,7 @@ md"""## Initialize Data Structures"""
 
 # ╔═╡ 66c95828-227c-4db5-a6f1-3e3004a99785
 begin
-	P,D=OCCA_FlowFields.setup(nmax=5)
+	P,D=OCCAmodule.setup(nmax=5)
 	"Done with FlowFields"
 end
 
@@ -43,7 +38,7 @@ end
 # ╔═╡ f199f321-976a-4ccd-a003-140211aa67fe
 begin	
 	I=Individuals(P,df.x,df.y,df.z,df.fid,
-		(🔴=OCCA_FlowFields.custom🔴,🔧=OCCA_FlowFields.custom🔧, D=D))
+		(🔴=OCCAmodule.custom🔴,🔧=OCCAmodule.custom🔧, D=D))
 	"Done with Individuals"
 end
 
diff --git a/ext/IndividualDisplacementsMITgcmExt.jl b/ext/IndividualDisplacementsMITgcmExt.jl
new file mode 100644
index 00000000..3bc5b874
--- /dev/null
+++ b/ext/IndividualDisplacementsMITgcmExt.jl
@@ -0,0 +1,12 @@
+
+module IndividualDisplacementsMITgcmExt
+
+using MITgcm, IndividualDisplacements
+import IndividualDisplacements: read_data_ECCO, NetCDF
+
+function read_data_ECCO(m0,v0,path,grid,k)
+    MITgcm.read_nctiles(joinpath(path,v0),v0,grid,I=(:,:,k,m0))
+end
+read_data_ECCO(m0,v0,P,D)=read_data_ECCO(m0,v0, joinpath(D.pth,v0) , P.u0.grid , D.k )
+
+end
diff --git a/src/IndividualDisplacements.jl b/src/IndividualDisplacements.jl
index e23de6e5..4ed0d999 100644
--- a/src/IndividualDisplacements.jl
+++ b/src/IndividualDisplacements.jl
@@ -1,7 +1,10 @@
 module IndividualDisplacements
 
 using MeshArrays, CyclicArrays, OrdinaryDiffEq, DataFrames, Random
+import NetCDF, CSV
+
 function data_path end
+function read_data_ECCO end
 
 include("API.jl")
 include("compute.jl")
@@ -9,6 +12,8 @@ include("data_wrangling.jl")
 include("toy_models.jl")
 include("initial_positions.jl")
 include("Downloads.jl")
+include("example_ECCO.jl")
+include("example_OCCA.jl")
 
 DiffEqBase.solve!(I::Individuals,args...)=∫!(I::Individuals,args...)
 DataFrames.groupby(I::Individuals,args...) = groupby(I.🔴,args...)
diff --git a/examples/worldwide/ECCO_FlowFields.jl b/src/example_ECCO.jl
similarity index 93%
rename from examples/worldwide/ECCO_FlowFields.jl
rename to src/example_ECCO.jl
index af4832a6..f29d85f6 100644
--- a/examples/worldwide/ECCO_FlowFields.jl
+++ b/src/example_ECCO.jl
@@ -1,11 +1,10 @@
-module ECCO_FlowFields
+module ECCO
 
-using IndividualDisplacements
 import MeshArrays, DataDeps, CSV, JLD2
 
 import IndividualDisplacements.OrdinaryDiffEq: solve, Tsit5, ODEProblem
-import IndividualDisplacements: update_location!
-import IndividualDisplacements: data_path
+import IndividualDisplacements: update_location!, Individuals, uvMeshArrays, uvwMeshArrays
+import IndividualDisplacements: FlowFields, data_path, read_data_ECCO
 import IndividualDisplacements.DataFrames: DataFrame
 import IndividualDisplacements.MeshArrays as MeshArrays
 import IndividualDisplacements.MeshArrays: gcmgrid, MeshArray, exchange
@@ -14,12 +13,6 @@ export init_FlowFields, init_storage
 export custom∫, custom🔧, custom🔴, custom∫!
 #export reset_📌!, init_z_if_needed
 
-import MITgcm, NetCDF
-function read_data(m0,v0,path,grid,k)
-    MITgcm.read_nctiles(joinpath(path,v0),v0,grid,I=(:,:,k,m0))
-end
-read_data(m0,v0,P,D)=read_data(m0,v0, joinpath(D.pth,v0) , P.u0.grid , D.k )
-
 """
     reset_📌!(I::Individuals,frac::Number,📌::Array)
 
@@ -150,19 +143,19 @@ function update_FlowFields!(P::uvMeshArrays,D::NamedTuple,t::AbstractFloat)
     P.v0[:]=Float32.(v0.MA[:])
     P.v1[:]=Float32.(v1.MA[:])
 
-    θ0=read_data(m0,"THETA",P,D)
+    θ0=read_data_ECCO(m0,"THETA",P,D)
     θ0[findall(isnan.(θ0))]=0.0 #mask with 0s rather than NaNs
     D.θ0[:]=Float32.(θ0[:,1])
 
-    θ1=read_data(m1,"THETA",P,D)
+    θ1=read_data_ECCO(m1,"THETA",P,D)
     θ1[findall(isnan.(θ1))]=0.0 #mask with 0s rather than NaNs
     D.θ1[:]=Float32.(θ1[:,1])
 
-    S0=read_data(m0,"SALT",P,D)
+    S0=read_data_ECCO(m0,"SALT",P,D)
     S0[findall(isnan.(S0))]=0.0 #mask with 0s rather than NaNs
     D.S0[:]=Float32.(S0[:,1])
 
-    S1=read_data(m1,"SALT",P,D)
+    S1=read_data_ECCO(m1,"SALT",P,D)
     S1[findall(isnan.(S1))]=0.0 #mask with 0s rather than NaNs
     D.S1[:]=Float32.(S1[:,1])
 
@@ -216,7 +209,7 @@ function update_FlowFields!(P::uvwMeshArrays,D::NamedTuple,t::AbstractFloat)
         u0[:,k]=tmpu.MA
         v0[:,k]=tmpv.MA
     end
-    w0=velocity_factor*read_data(m0,"WVELMASS",joinpath(D.pth,"WVELMASS"),P.u0.grid,:)
+    w0=velocity_factor*read_data_ECCO(m0,"WVELMASS",joinpath(D.pth,"WVELMASS"),P.u0.grid,:)
     w0[findall(isnan.(w0))]=0.0 #mask with 0s rather than NaNs
 
     (U,V)=read_velocities(P.u0.grid,m1,D.pth)
@@ -228,7 +221,7 @@ function update_FlowFields!(P::uvwMeshArrays,D::NamedTuple,t::AbstractFloat)
         u1[:,k]=tmpu.MA
         v1[:,k]=tmpv.MA
     end
-    w1=velocity_factor*read_data(m1,"WVELMASS",joinpath(D.pth,"WVELMASS"),P.u0.grid,:)
+    w1=velocity_factor*read_data_ECCO(m1,"WVELMASS",joinpath(D.pth,"WVELMASS"),P.u0.grid,:)
     w1[findall(isnan.(w1))]=0.0 #mask with 0s rather than NaNs
 
     P.u0[:,:]=Float32.(u0[:,:])
@@ -246,19 +239,19 @@ function update_FlowFields!(P::uvwMeshArrays,D::NamedTuple,t::AbstractFloat)
     P.w0[:,nr+1]=0*exchange(-w0[:,1],1).MA
     P.w1[:,nr+1]=0*exchange(-w1[:,1],1).MA
 
-    θ0=read_data(m0,"THETA",joinpath(D.pth,"THETA"),P.u0.grid,:)
+    θ0=read_data_ECCO(m0,"THETA",joinpath(D.pth,"THETA"),P.u0.grid,:)
     θ0[findall(isnan.(θ0))]=0.0 #mask with 0s rather than NaNs
     D.θ0[:,:]=Float32.(θ0[:,:])
 
-    θ1=read_data(m1,"THETA",joinpath(D.pth,"THETA"),P.u0.grid,:)
+    θ1=read_data_ECCO(m1,"THETA",joinpath(D.pth,"THETA"),P.u0.grid,:)
     θ1[findall(isnan.(θ1))]=0.0 #mask with 0s rather than NaNs
     D.θ1[:,:]=Float32.(θ1[:,:])
 
-    S0=read_data(m0,"SALT",joinpath(D.pth,"SALT"),P.u0.grid,:)
+    S0=read_data_ECCO(m0,"SALT",joinpath(D.pth,"SALT"),P.u0.grid,:)
     S0[findall(isnan.(S0))]=0.0 #mask with 0s rather than NaNs
     D.S0[:,:]=Float32.(S0[:,:])
 
-    S1=read_data(m1,"SALT",joinpath(D.pth,"SALT"),P.u0.grid,:)
+    S1=read_data_ECCO(m1,"SALT",joinpath(D.pth,"SALT"),P.u0.grid,:)
     S1[findall(isnan.(S1))]=0.0 #mask with 0s rather than NaNs
     D.S1[:,:]=Float32.(S1[:,:])
 
@@ -293,8 +286,8 @@ end
 Read velocity components `u,v` from files in `pth`for time `t`
 """
 function read_velocities(γ::gcmgrid,t::Int,pth::String)
-    u=read_data(t,"UVELMASS",joinpath(pth,"UVELMASS"),γ,:)
-    v=read_data(t,"VVELMASS",joinpath(pth,"VVELMASS"),γ,:)
+    u=read_data_ECCO(t,"UVELMASS",joinpath(pth,"UVELMASS"),γ,:)
+    v=read_data_ECCO(t,"VVELMASS",joinpath(pth,"VVELMASS"),γ,:)
     return u,v
 end
 
diff --git a/examples/worldwide/OCCA_FlowFields.jl b/src/example_OCCA.jl
similarity index 96%
rename from examples/worldwide/OCCA_FlowFields.jl
rename to src/example_OCCA.jl
index 9d4c623e..a5e99a62 100644
--- a/examples/worldwide/OCCA_FlowFields.jl
+++ b/src/example_OCCA.jl
@@ -1,7 +1,7 @@
-module OCCA_FlowFields
+module OCCA
 
-using IndividualDisplacements
-import IndividualDisplacements: data_path
+import IndividualDisplacements: data_path, uvwMeshArrays, FlowFields
+import IndividualDisplacements: postprocess_MeshArray, add_lonlat!, interp_to_xy
 
 ##
 
@@ -168,4 +168,4 @@ function custom∫(prob)
          trajectories=length(u0),saveat=365/12*86400.0)
 end
 
-end #module OCCA_FlowFields
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index fe348036..8b43d239 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,5 +1,6 @@
 using Test, Documenter
 using IndividualDisplacements, Climatology, MeshArrays, NetCDF, Suppressor
+import MITgcm
 
 Climatology.get_ecco_velocity_if_needed()
 Climatology.get_occa_velocity_if_needed()
@@ -17,11 +18,10 @@ end
 
 @testset "global" begin
     p0=IndividualDisplacements.datadeps.getdata("global_ocean_circulation_inputs")
-    p1=dirname(pathof(IndividualDisplacements))
-    include(joinpath(p1,"../examples/worldwide/ECCO_FlowFields.jl"))
-    P,D=ECCO_FlowFields.init_FlowFields()
+    ECCOmodule=IndividualDisplacements.ECCO
+    P,D=ECCOmodule.init_FlowFields()
     file_input=joinpath(p0,"initial_10_1.csv")
-    df = ECCO_FlowFields.init_positions(10,filename=file_input)
+    df = IndividualDisplacements.init.init_positions(10,filename=file_input)
     I=Individuals(P,df.x,df.y,df.f,(D=D,))
     T=(0.0,I.P.T[2])
     ∫!(I,T)