diff --git a/src/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl b/src/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl
index 29cd770..640e2d2 100644
--- a/src/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl
+++ b/src/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl
@@ -1,9 +1,9 @@
 function JetProp2DAcoTTIDenQ_DEO2_FDTD(;
-        v = Float32[],
-        b = Float32[],
-        ϵ = Float32[],
-        η = Float32[],
-        θ = Float32[],
+        v::Array{Float32,2} = ones(Float32,0,0),
+        b::Array{Float32,2} = ones(Float32,0,0),
+        ϵ::Array{Float32,2} = ones(Float32,0,0),
+        η::Array{Float32,2} = ones(Float32,0,0),
+        θ::Array{Float32,2} = ones(Float32,0,0),
         f = 0.85,
         srcfieldfile = joinpath(tempdir(), "field-$(uuid4()).bin"),
         comptype = nothing,
@@ -39,6 +39,7 @@ function JetProp2DAcoTTIDenQ_DEO2_FDTD(;
         imgcondition = "standard",
         nthreads = Sys.CPU_THREADS,
         reportinterval = 500)
+        
     # active and passive earth model properties.  The active set is in the model-space
     active_modelset = Dict{String,Int}()
     passive_modelset = Dict{String,Array{Float32,2}}()
@@ -106,7 +107,7 @@ function JetProp2DAcoTTIDenQ_DEO2_FDTD(;
     C = WaveFD.comptype(comptype, Float32)[1]
     compressor = Dict{String,WaveFD.Compressor{Float32,Float32,C,2}}()
 
-    # we need to serialize for the data, but not for the linearization
+    # we need to serialize "pold" for the data but, not (necessarily) for the linearization
     _active_wavefields = "pold" ∈ active_wavefields ? active_wavefields : [active_wavefields;"pold"]
 
     for active_wavefield in _active_wavefields
@@ -124,7 +125,6 @@ function JetProp2DAcoTTIDenQ_DEO2_FDTD(;
         error("Supplied imaging condition 'imgcondition' is not in [standard, FWI, RTM]")
     end
         
-    # construct:
     Jet(
         dom = dom,
         rng = rng,
@@ -441,10 +441,10 @@ function JopProp2DAcoTTIDenQ_DEO2_FDTD_nonlinearforward!(d::AbstractArray, m::Ab
         qMin = kwargs[:qMin],
         qInterior = kwargs[:qInterior])
 
-    # wave-fields
+    # wavefields
     wavefields = Dict("pcur"=>WaveFD.PCur(p), "pold"=>WaveFD.POld(p), "pspace"=>WaveFD.PSpace(p), "mcur"=>WaveFD.MCur(p), "mold"=>WaveFD.MOld(p), "mspace"=>WaveFD.MSpace(p))
 
-    # earth model
+    # ginsu'd earth model
     model_ginsu = Dict("v"=>WaveFD.V(p), "b"=>WaveFD.B(p), "ϵ"=>WaveFD.Eps(p), "η"=>WaveFD.Eta(p), "f"=>WaveFD.F(p), 
         "cosθ"=>WaveFD.CosTheta(p), "sinθ"=>WaveFD.SinTheta(p))
 
@@ -452,10 +452,14 @@ function JopProp2DAcoTTIDenQ_DEO2_FDTD_nonlinearforward!(d::AbstractArray, m::Ab
     for prop in keys(kwargs[:active_modelset])
         sub!(model_ginsu[prop], kwargs[:ginsu], @view(m[:,:,kwargs[:active_modelset][prop]]), extend=true)
     end
+    # ginsu'd earth model (passive-set)
     for prop in keys(kwargs[:passive_modelset])
         sub!(model_ginsu[prop], kwargs[:ginsu], kwargs[:passive_modelset][prop], extend=true)
     end
+
+    # TODO why no ginsu for f?
     model_ginsu["f"] .= kwargs[:f]
+
     ginsu_interior_range = interior(kwargs[:ginsu])
 
     # we need to serialize "pold" for the data but, not (necessarily) for the linearization
@@ -500,10 +504,13 @@ function JopProp2DAcoTTIDenQ_DEO2_FDTD_nonlinearforward!(d::AbstractArray, m::Ab
     if kwargs[:srcfieldfile] != ""
         for active_wavefield in active_wavefields
             filename = "$(kwargs[:srcfieldfile])-$(active_wavefield)"
-            if isfile(filename) == true
-                rm(filename)
+            try
+                isfile(filename) && rm(filename)
+                iofield[active_wavefield] = open(filename, "w")
+            catch
+                @info "Unable to open $(filename) on proc $(myid()) - $(gethostname())"
+                rethrow()
             end
-            iofield[active_wavefield] = open(filename, "w")
             open(kwargs[:compressor][active_wavefield])
         end
     end
@@ -516,7 +523,8 @@ function JopProp2DAcoTTIDenQ_DEO2_FDTD_nonlinearforward!(d::AbstractArray, m::Ab
     set_zero_subnormals(true)
     for it = 1:ntmod_wav
         if kwargs[:reportinterval] != 0 && (it % kwargs[:reportinterval] == 0 || it == ntmod_wav)
-            JopProp2DAcoTTIDenQ_DEO2_FDTD_write_history_nl(kwargs[:ginsu], it, ntmod_wav, time()-time1, cumtime_io, cumtime_ex, cumtime_pr, wavefields["pcur"], d)
+            JopProp2DAcoTTIDenQ_DEO2_FDTD_write_history_nl(kwargs[:ginsu], it, ntmod_wav, time()-time1, 
+                cumtime_io, cumtime_ex, cumtime_pr, wavefields["pcur"], d)
         end
 
         # propagate and wavefield swap
@@ -542,9 +550,11 @@ function JopProp2DAcoTTIDenQ_DEO2_FDTD_nonlinearforward!(d::AbstractArray, m::Ab
             if kwargs[:srcfieldfile] != ""
                 cumtime_io += @elapsed for active_wavefield in active_wavefields
                     if kwargs[:isinterior]
-                        WaveFD.compressedwrite(iofield[active_wavefield],  kwargs[:compressor][active_wavefield], div(it-1,itskip)+1, wavefields[active_wavefield], ginsu_interior_range)
+                        WaveFD.compressedwrite(iofield[active_wavefield], kwargs[:compressor][active_wavefield],
+                            div(it-1,itskip)+1, wavefields[active_wavefield], ginsu_interior_range)
                     else
-                        WaveFD.compressedwrite(iofield[active_wavefield],  kwargs[:compressor][active_wavefield], div(it-1,itskip)+1, wavefields[active_wavefield])
+                        WaveFD.compressedwrite(iofield[active_wavefield], kwargs[:compressor][active_wavefield], 
+                            div(it-1,itskip)+1, wavefields[active_wavefield])
                     end
                 end
             end
@@ -574,15 +584,15 @@ function JopProp2DAcoTTIDenQ_DEO2_FDTD_f!(d::AbstractArray, m::AbstractArray{Flo
         kwargs[:srcfieldhost][] = gethostname()
     else
         field = Array{Float32}(undef,size(kwargs[:ginsu], interior=kwargs[:isinterior]))
-        local iz, ix, ic
+        local iz, ix, c
         if kwargs[:interpmethod] == :hicks
             iz, ix, c = WaveFD.hickscoeffs(kwargs[:dz], kwargs[:dx], 
-                origin(kwargs[:ginsu],interior=kwargs[:isinterior])..., 
-                size(kwargs[:ginsu],interior=kwargs[:isinterior])..., kwargs[:rz], kwargs[:rx])
+                origin(kwargs[:ginsu], interior=kwargs[:isinterior])..., 
+                size(kwargs[:ginsu], interior=kwargs[:isinterior])..., kwargs[:rz], kwargs[:rx])
         else
             iz, ix, c = WaveFD.linearcoeffs(kwargs[:dz], kwargs[:dx], 
-                origin(kwargs[:ginsu],interior=kwargs[:isinterior])..., 
-                size(kwargs[:ginsu],interior=kwargs[:isinterior])..., kwargs[:rz], kwargs[:rx])
+                origin(kwargs[:ginsu], interior=kwargs[:isinterior])..., 
+                size(kwargs[:ginsu], interior=kwargs[:isinterior])..., kwargs[:rz], kwargs[:rx])
         end
 
         iofield = open("$(kwargs[:srcfieldfile])-pold")
@@ -598,8 +608,8 @@ function JopProp2DAcoTTIDenQ_DEO2_FDTD_f!(d::AbstractArray, m::AbstractArray{Flo
 end
 
 function JopProp2DAcoTTIDenQ_DEO2_FDTD_df!(δd::AbstractArray, δm::AbstractArray; kwargs...)
-    nz_ginsu,nx_ginsu=size(kwargs[:ginsu])
-    z0_ginsu,x0_ginsu=origin(kwargs[:ginsu])
+    nz_ginsu,nx_ginsu = size(kwargs[:ginsu])
+    z0_ginsu,x0_ginsu = origin(kwargs[:ginsu])
     p = Prop2DAcoTTIDenQ_DEO2_FDTD(
         nz = nz_ginsu,
         nx = nx_ginsu,
@@ -615,21 +625,22 @@ function JopProp2DAcoTTIDenQ_DEO2_FDTD_df!(δd::AbstractArray, δm::AbstractArra
         qMin = kwargs[:qMin],
         qInterior = kwargs[:qInterior])
 
-    # wave-fields
+    # wavefields
     pcur,pold = WaveFD.PCur(p),WaveFD.POld(p)
 
-    # earth model
-    earth_ginsu = Dict("v"=>WaveFD.V(p), "b"=>WaveFD.B(p), "ϵ"=>WaveFD.Eps(p), "η"=>WaveFD.Eta(p), "f"=>WaveFD.F(p), 
+    # ginsu'd earth model
+    model_ginsu = Dict("v"=>WaveFD.V(p), "b"=>WaveFD.B(p), "ϵ"=>WaveFD.Eps(p), "η"=>WaveFD.Eta(p), "f"=>WaveFD.F(p), 
         "cosθ"=>WaveFD.CosTheta(p), "sinθ"=>WaveFD.SinTheta(p))
 
     # ginsu'd earth model (active-set)
     for prop in keys(kwargs[:active_modelset])
-        sub!(earth_ginsu[prop], kwargs[:ginsu], @view(kwargs[:mₒ][:,:,kwargs[:active_modelset][prop]]), extend=true)
+        sub!(model_ginsu[prop], kwargs[:ginsu], @view(kwargs[:mₒ][:,:,kwargs[:active_modelset][prop]]), extend=true)
     end
+    # ginsu'd earth model (passive-set)
     for prop in keys(kwargs[:passive_modelset])
-        sub!(earth_ginsu[prop], kwargs[:ginsu], kwargs[:passive_modelset][prop], extend=true)
+        sub!(model_ginsu[prop], kwargs[:ginsu], kwargs[:passive_modelset][prop], extend=true)
     end
-    earth_ginsu["f"] .= kwargs[:f]
+    model_ginsu["f"] .= kwargs[:f]
 
     δm_ginsu = Dict{String,Array{Float32,2}}()
     for prop in keys(kwargs[:active_modelset])
@@ -639,7 +650,7 @@ function JopProp2DAcoTTIDenQ_DEO2_FDTD_df!(δd::AbstractArray, δm::AbstractArra
     ginsu_interior_range = interior(kwargs[:ginsu])
 
     # pre-compute receiver interpolation coefficients
-    local iz,ix,c
+    local iz, ix, c
     if kwargs[:interpmethod] == :hicks
         iz, ix, c = WaveFD.hickscoeffs(kwargs[:dz], kwargs[:dx], z0_ginsu, x0_ginsu, nz_ginsu, nx_ginsu, kwargs[:rz], kwargs[:rx])
     else
@@ -687,9 +698,11 @@ function JopProp2DAcoTTIDenQ_DEO2_FDTD_df!(δd::AbstractArray, δm::AbstractArra
             # read source field from disk
             cumtime_io += @elapsed for active_wavefield in kwargs[:active_wavefields]
                 if kwargs[:isinterior]
-                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield], div(it-1,itskip)+1, wavefields[active_wavefield], ginsu_interior_range)
+                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield], 
+                        div(it-1,itskip)+1, wavefields[active_wavefield], ginsu_interior_range)
                 else
-                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield], div(it-1,itskip)+1, wavefields[active_wavefield])
+                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield], 
+                        div(it-1,itskip)+1, wavefields[active_wavefield])
                 end
             end
 
@@ -735,22 +748,22 @@ function JopProp2DAcoTTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
         qMin = kwargs[:qMin],
         qInterior = kwargs[:qInterior])
 
-    # wave-fields
+    # wavefields
     pcur,pold = WaveFD.PCur(p),WaveFD.POld(p)
 
     # ginsu'd earth model
-    earth_ginsu = Dict("v"=>WaveFD.V(p), "b"=>WaveFD.B(p), "ϵ"=>WaveFD.Eps(p), "η"=>WaveFD.Eta(p), "f"=>WaveFD.F(p), 
+    model_ginsu = Dict("v"=>WaveFD.V(p), "b"=>WaveFD.B(p), "ϵ"=>WaveFD.Eps(p), "η"=>WaveFD.Eta(p), "f"=>WaveFD.F(p), 
         "cosθ"=>WaveFD.CosTheta(p), "sinθ"=>WaveFD.SinTheta(p))
 
-    # active model-set
+    # ginsu'd earth model (active-set)
     for prop in keys(kwargs[:active_modelset])
-        sub!(earth_ginsu[prop], kwargs[:ginsu], @view(kwargs[:mₒ][:,:,kwargs[:active_modelset][prop]]), extend=true)
+        sub!(model_ginsu[prop], kwargs[:ginsu], @view(kwargs[:mₒ][:,:,kwargs[:active_modelset][prop]]), extend=true)
     end
-    # passive model-set
+    # ginsu'd earth model (passive-set)
     for prop in keys(kwargs[:passive_modelset])
-        sub!(earth_ginsu[prop], kwargs[:ginsu], kwargs[:passive_modelset][prop], extend=true)
+        sub!(model_ginsu[prop], kwargs[:ginsu], kwargs[:passive_modelset][prop], extend=true)
     end
-    earth_ginsu["f"] .= kwargs[:f]
+    model_ginsu["f"] .= kwargs[:f]
 
     δm_ginsu = Dict{String,Array{Float32,2}}()
     for prop in keys(kwargs[:active_modelset])
@@ -760,7 +773,7 @@ function JopProp2DAcoTTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
     ginsu_interior_range = interior(kwargs[:ginsu])
 
     # Get receiver interpolation coefficients
-    local iz,ix,c
+    local iz, ix, c
     if kwargs[:interpmethod] == :hicks
         iz, ix, c = WaveFD.hickscoeffs(kwargs[:dz], kwargs[:dx], z0_ginsu, x0_ginsu, nz_ginsu, nx_ginsu, kwargs[:rz], kwargs[:rx])
     else
@@ -770,7 +783,7 @@ function JopProp2DAcoTTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
         for jx = 1:size(c[i], 2), jz = 1:size(c[i], 1)
             kz = iz[i][jz]
             kx = ix[i][jx]
-            c[i][jz,jx] *= kwargs[:dtmod]^2 * earth_ginsu["v"][kz,kx]^2 / earth_ginsu["b"][kz,kx]
+            c[i][jz,jx] *= kwargs[:dtmod]^2 * model_ginsu["v"][kz,kx]^2 / model_ginsu["b"][kz,kx]
         end
     end
     blks = WaveFD.source_blocking(nz_ginsu, nx_ginsu, kwargs[:nbz_inject], kwargs[:nbx_inject], iz, ix, c)
@@ -806,7 +819,8 @@ function JopProp2DAcoTTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
     set_zero_subnormals(true)
     for it = ntmod:-1:1
         if kwargs[:reportinterval] != 0 && (it % kwargs[:reportinterval] == 0 || it == ntmod)
-            JopProp2DAcoTTIDenQ_DEO2_FDTD_write_history_ln(kwargs[:ginsu], it, ntmod, time()-time1, cumtime_io, cumtime_ex, cumtime_im, cumtime_pr, pcur, δdinterp, "adjoint")
+            JopProp2DAcoTTIDenQ_DEO2_FDTD_write_history_ln(kwargs[:ginsu], it, ntmod, time()-time1,
+                cumtime_io, cumtime_ex, cumtime_im, cumtime_pr, pcur, δdinterp, "adjoint")
         end
 
         # propagate and wavefield swap
@@ -820,9 +834,11 @@ function JopProp2DAcoTTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
             # read source field from disk
             cumtime_io += @elapsed for active_wavefield in kwargs[:active_wavefields]
                 if kwargs[:isinterior]
-                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield], div(it-1,itskip)+1, wavefields[active_wavefield], ginsu_interior_range)
+                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield],
+                        div(it-1,itskip)+1, wavefields[active_wavefield], ginsu_interior_range)
                 else
-                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield], div(it-1,itskip)+1, wavefields[active_wavefield])
+                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield],
+                        div(it-1,itskip)+1, wavefields[active_wavefield])
                 end
             end
 
@@ -831,7 +847,8 @@ function JopProp2DAcoTTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
         end
     end
     set_zero_subnormals(false)
-    JopProp2DAcoTTIDenQ_DEO2_FDTD_stats(kwargs[:stats], kwargs[:ginsu], ntmod, time()-time1, cumtime_io, cumtime_ex, cumtime_im)
+    JopProp2DAcoTTIDenQ_DEO2_FDTD_stats(kwargs[:stats], kwargs[:ginsu], ntmod, time()-time1,
+        cumtime_io, cumtime_ex, cumtime_im)
 
     # undo ginsu
     for prop in keys(kwargs[:active_modelset])
@@ -848,11 +865,11 @@ function JopProp2DAcoTTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
     δm
 end
 
-modelindex(F::Jop{T}, key::AbstractString) where {D,R,T<:Jet{D,R,typeof(JopProp2DAcoTTIDenQ_DEO2_FDTD_f!)}} = state(F).active_modelset[key]
+modelindex(F::Jop{T}, key) where {D,R,T<:Jet{D,R,typeof(JopProp2DAcoTTIDenQ_DEO2_FDTD_f!)}} = state(F).active_modelset[key]
 
-function srcillum!(γ, A::T, m::AbstractArray{Float32}) where {D,R,J<:Jet{D,R,typeof(JopProp2DAcoTTIDenQ_DEO2_FDTD_f!)},T<:Jop{J}}
+function srcillum!(γ, A::Jop{T}, m::AbstractArray{Float32}) where {D,R,T<:Jet{D,R,typeof(JopProp2DAcoTTIDenQ_DEO2_FDTD_f!)}}
     s = state(A)
-    isvalid, _chksum = isvalid_srcfieldfile(m, s.srcfieldhost[], s.srcfieldfile*"-pold", s.chksum[])
+    isvalid, _chksum = isvalid_srcfieldfile(jet(A).mₒ, s.srcfieldhost[], s.srcfieldfile*"-pold", s.chksum[])
     if !isvalid
         JopProp2DAcoTTIDenQ_DEO2_FDTD_nonlinearforward!(Array{Float32}(undef,0,0), m; s...)
         s.chksum[] = _chksum
@@ -903,7 +920,7 @@ end
 
 function Base.close(jet::Jet{D,R,typeof(JopProp2DAcoTTIDenQ_DEO2_FDTD_f!)}) where {D,R}
     rm("$(state(jet).srcfieldfile)-pold", force=true)
-    rm("$(state(jet).srcfieldfile)-mold", force=true)
     rm("$(state(jet).srcfieldfile)-pspace", force=true)
+    rm("$(state(jet).srcfieldfile)-mold", force=true)
     rm("$(state(jet).srcfieldfile)-mspace", force=true)
 end
diff --git a/src/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl b/src/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl
index 2d394c1..acb44c6 100644
--- a/src/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl
+++ b/src/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl
@@ -1,8 +1,8 @@
 function JetProp2DAcoVTIDenQ_DEO2_FDTD(;
-        v = Float32[],
-        b = Float32[],
-        ϵ = Float32[],
-        η = Float32[],
+        v::Array{Float32,2} = ones(Float32,0,0),
+        b::Array{Float32,2} = ones(Float32,0,0),
+        ϵ::Array{Float32,2} = ones(Float32,0,0),
+        η::Array{Float32,2} = ones(Float32,0,0),
         f = 0.85,
         srcfieldfile = joinpath(tempdir(), "field-$(uuid4()).bin"),
         comptype = nothing,
@@ -38,6 +38,7 @@ function JetProp2DAcoVTIDenQ_DEO2_FDTD(;
         imgcondition = "standard",
         nthreads = Sys.CPU_THREADS,
         reportinterval = 500)
+
     # active and passive earth model properties.  The active set is in the model-space
     active_modelset = Dict{String,Int}()
     passive_modelset = Dict{String,Array{Float32,2}}()
@@ -103,7 +104,7 @@ function JetProp2DAcoVTIDenQ_DEO2_FDTD(;
     C = WaveFD.comptype(comptype, Float32)[1]
     compressor = Dict{String,WaveFD.Compressor{Float32,Float32,C,2}}()
 
-    # we need to serialize for the data, but not for the linearization
+    # we need to serialize "pold" for the data but, not (necessarily) for the linearization
     _active_wavefields = "pold" ∈ active_wavefields ? active_wavefields : [active_wavefields;"pold"]
 
     for active_wavefield in _active_wavefields
@@ -433,19 +434,22 @@ function JopProp2DAcoVTIDenQ_DEO2_FDTD_nonlinearforward!(d::AbstractArray, m::Ab
         qMin = kwargs[:qMin],
         qInterior = kwargs[:qInterior])
 
-    # wave-fields
+    # wavefields
     wavefields = Dict("pcur"=>WaveFD.PCur(p), "pold"=>WaveFD.POld(p), "pspace"=>WaveFD.PSpace(p), "mcur"=>WaveFD.MCur(p), "mold"=>WaveFD.MOld(p), "mspace"=>WaveFD.MSpace(p))
 
-    # earth model
+    # ginsu'd earth model
     model_ginsu = Dict("v"=>WaveFD.V(p), "b"=>WaveFD.B(p), "ϵ"=>WaveFD.Eps(p), "η"=>WaveFD.Eta(p), "f"=>WaveFD.F(p))
 
     # ginsu'd earth model (active-set)
     for prop in keys(kwargs[:active_modelset])
         sub!(model_ginsu[prop], kwargs[:ginsu], @view(m[:,:,kwargs[:active_modelset][prop]]), extend=true)
     end
+    # ginsu'd earth model (passive-set)
     for prop in keys(kwargs[:passive_modelset])
         sub!(model_ginsu[prop], kwargs[:ginsu], kwargs[:passive_modelset][prop], extend=true)
     end
+
+    # TODO why no ginsu for f?
     model_ginsu["f"] .= kwargs[:f]
 
     ginsu_interior_range = interior(kwargs[:ginsu])
@@ -493,9 +497,7 @@ function JopProp2DAcoVTIDenQ_DEO2_FDTD_nonlinearforward!(d::AbstractArray, m::Ab
         for active_wavefield in active_wavefields
             filename = "$(kwargs[:srcfieldfile])-$(active_wavefield)"
             try
-                if isfile(filename) == true
-                    rm(filename)
-                end
+                isfile(filename) && rm(filename)
                 iofield[active_wavefield] = open(filename, "w")
             catch
                 @info "Unable to open $(filename) on proc $(myid()) - $(gethostname())"
@@ -513,7 +515,8 @@ function JopProp2DAcoVTIDenQ_DEO2_FDTD_nonlinearforward!(d::AbstractArray, m::Ab
     set_zero_subnormals(true)
     for it = 1:ntmod_wav
         if kwargs[:reportinterval] != 0 && (it % kwargs[:reportinterval] == 0 || it == ntmod_wav)
-            JopProp2DAcoVTIDenQ_DEO2_FDTD_write_history_nl(kwargs[:ginsu], it, ntmod_wav, time()-time1, cumtime_io, cumtime_ex, cumtime_pr, wavefields["pcur"], d)
+            JopProp2DAcoVTIDenQ_DEO2_FDTD_write_history_nl(kwargs[:ginsu], it, ntmod_wav, time()-time1, 
+                cumtime_io, cumtime_ex, cumtime_pr, wavefields["pcur"], d)
         end
 
         # propagate and wavefield swap
@@ -539,11 +542,11 @@ function JopProp2DAcoVTIDenQ_DEO2_FDTD_nonlinearforward!(d::AbstractArray, m::Ab
             if kwargs[:srcfieldfile] != ""
                 cumtime_io += @elapsed for active_wavefield in active_wavefields
                     if kwargs[:isinterior]
-                        WaveFD.compressedwrite(iofield[active_wavefield], kwargs[:compressor][active_wavefield], div(it-1,itskip)+1,
-                            wavefields[active_wavefield], ginsu_interior_range)
+                        WaveFD.compressedwrite(iofield[active_wavefield], kwargs[:compressor][active_wavefield], 
+                            div(it-1,itskip)+1, wavefields[active_wavefield], ginsu_interior_range)
                     else
-                        WaveFD.compressedwrite(iofield[active_wavefield], kwargs[:compressor][active_wavefield], div(it-1,itskip)+1,
-                            wavefields[active_wavefield])
+                        WaveFD.compressedwrite(iofield[active_wavefield], kwargs[:compressor][active_wavefield], 
+                            div(it-1,itskip)+1, wavefields[active_wavefield])
                     end
                 end
             end
@@ -614,20 +617,21 @@ function JopProp2DAcoVTIDenQ_DEO2_FDTD_df!(δd::AbstractArray, δm::AbstractArra
         qMin = kwargs[:qMin],
         qInterior = kwargs[:qInterior])
 
-    # wave-fields
+    # wavefields
     pcur,pold = WaveFD.PCur(p),WaveFD.POld(p)
 
-    # earth model
-    earth_ginsu = Dict("v"=>WaveFD.V(p), "b"=>WaveFD.B(p), "ϵ"=>WaveFD.Eps(p), "η"=>WaveFD.Eta(p), "f"=>WaveFD.F(p))
+    # ginsu'd earth model
+    model_ginsu = Dict("v"=>WaveFD.V(p), "b"=>WaveFD.B(p), "ϵ"=>WaveFD.Eps(p), "η"=>WaveFD.Eta(p), "f"=>WaveFD.F(p))
 
     # ginsu'd earth model (active-set)
     for prop in keys(kwargs[:active_modelset])
-        sub!(earth_ginsu[prop], kwargs[:ginsu], @view(kwargs[:mₒ][:,:,kwargs[:active_modelset][prop]]), extend=true)
+        sub!(model_ginsu[prop], kwargs[:ginsu], @view(kwargs[:mₒ][:,:,kwargs[:active_modelset][prop]]), extend=true)
     end
+    # ginsu'd earth model (passive-set)
     for prop in keys(kwargs[:passive_modelset])
-        sub!(earth_ginsu[prop], kwargs[:ginsu], kwargs[:passive_modelset][prop], extend=true)
+        sub!(model_ginsu[prop], kwargs[:ginsu], kwargs[:passive_modelset][prop], extend=true)
     end
-    earth_ginsu["f"] .= kwargs[:f]
+    model_ginsu["f"] .= kwargs[:f]
 
     δm_ginsu = Dict{String,Array{Float32,2}}()
     for prop in keys(kwargs[:active_modelset])
@@ -637,7 +641,7 @@ function JopProp2DAcoVTIDenQ_DEO2_FDTD_df!(δd::AbstractArray, δm::AbstractArra
     ginsu_interior_range = interior(kwargs[:ginsu])
 
     # pre-compute receiver interpolation coefficients
-    local iz,ix,c
+    local iz, ix, c
     if kwargs[:interpmethod] == :hicks
         iz, ix, c = WaveFD.hickscoeffs(kwargs[:dz], kwargs[:dx], z0_ginsu, x0_ginsu, nz_ginsu, nx_ginsu, kwargs[:rz], kwargs[:rx])
     else
@@ -685,11 +689,11 @@ function JopProp2DAcoVTIDenQ_DEO2_FDTD_df!(δd::AbstractArray, δm::AbstractArra
             # read source field from disk
             cumtime_io += @elapsed for active_wavefield in kwargs[:active_wavefields]
                 if kwargs[:isinterior]
-                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield], div(it-1,itskip)+1,
-                        wavefields[active_wavefield], ginsu_interior_range)
+                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield],
+                        div(it-1,itskip)+1, wavefields[active_wavefield], ginsu_interior_range)
                 else
-                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield], div(it-1,itskip)+1,
-                        wavefields[active_wavefield])
+                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield],
+                        div(it-1,itskip)+1, wavefields[active_wavefield])
                 end
             end
 
@@ -735,21 +739,21 @@ function JopProp2DAcoVTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
         qMin = kwargs[:qMin],
         qInterior = kwargs[:qInterior])
 
-    # wave-fields
+    # wavefields
     pcur,pold = WaveFD.PCur(p),WaveFD.POld(p)
 
     # ginsu'd earth model
-    earth_ginsu = Dict("v"=>WaveFD.V(p), "b"=>WaveFD.B(p), "ϵ"=>WaveFD.Eps(p), "η"=>WaveFD.Eta(p), "f"=>WaveFD.F(p))
+    model_ginsu = Dict("v"=>WaveFD.V(p), "b"=>WaveFD.B(p), "ϵ"=>WaveFD.Eps(p), "η"=>WaveFD.Eta(p), "f"=>WaveFD.F(p))
 
-    # active model-set
+    # ginsu'd earth model (active-set)
     for prop in keys(kwargs[:active_modelset])
-        sub!(earth_ginsu[prop], kwargs[:ginsu], @view(kwargs[:mₒ][:,:,kwargs[:active_modelset][prop]]), extend=true)
+        sub!(model_ginsu[prop], kwargs[:ginsu], @view(kwargs[:mₒ][:,:,kwargs[:active_modelset][prop]]), extend=true)
     end
-    # passive model-set
+    # ginsu'd earth model (passive-set)
     for prop in keys(kwargs[:passive_modelset])
-        sub!(earth_ginsu[prop], kwargs[:ginsu], kwargs[:passive_modelset][prop], extend=true)
+        sub!(model_ginsu[prop], kwargs[:ginsu], kwargs[:passive_modelset][prop], extend=true)
     end
-    earth_ginsu["f"] .= kwargs[:f]
+    model_ginsu["f"] .= kwargs[:f]
 
     δm_ginsu = Dict{String,Array{Float32,2}}()
     for prop in keys(kwargs[:active_modelset])
@@ -769,7 +773,7 @@ function JopProp2DAcoVTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
         for jx = 1:size(c[i], 2), jz = 1:size(c[i], 1)
             kz = iz[i][jz]
             kx = ix[i][jx]
-            c[i][jz,jx] *= kwargs[:dtmod]^2 * earth_ginsu["v"][kz,kx]^2 / earth_ginsu["b"][kz,kx]
+            c[i][jz,jx] *= kwargs[:dtmod]^2 * model_ginsu["v"][kz,kx]^2 / model_ginsu["b"][kz,kx]
         end
     end
     blks = WaveFD.source_blocking(nz_ginsu, nx_ginsu, kwargs[:nbz_inject], kwargs[:nbx_inject], iz, ix, c)
@@ -805,7 +809,8 @@ function JopProp2DAcoVTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
     set_zero_subnormals(true)
     for it = ntmod:-1:1
         if kwargs[:reportinterval] != 0 && (it % kwargs[:reportinterval] == 0 || it == ntmod)
-            JopProp2DAcoVTIDenQ_DEO2_FDTD_write_history_ln(kwargs[:ginsu], it, ntmod, time()-time1, cumtime_io, cumtime_ex, cumtime_im, cumtime_pr, pcur, δdinterp, "adjoint")
+            JopProp2DAcoVTIDenQ_DEO2_FDTD_write_history_ln(kwargs[:ginsu], it, ntmod, time()-time1,
+                cumtime_io, cumtime_ex, cumtime_im, cumtime_pr, pcur, δdinterp, "adjoint")
         end
 
         # propagate and wavefield swap
@@ -819,11 +824,11 @@ function JopProp2DAcoVTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
             # read source field from disk
             cumtime_io += @elapsed for active_wavefield in kwargs[:active_wavefields]
                 if kwargs[:isinterior]
-                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield], div(it-1,itskip)+1,
-                        wavefields[active_wavefield], ginsu_interior_range)
+                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield],
+                        div(it-1,itskip)+1, wavefields[active_wavefield], ginsu_interior_range)
                 else
-                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield], div(it-1,itskip)+1,
-                        wavefields[active_wavefield])
+                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield],
+                        div(it-1,itskip)+1, wavefields[active_wavefield])
                 end
             end
 
@@ -832,7 +837,8 @@ function JopProp2DAcoVTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
         end
     end
     set_zero_subnormals(false)
-    JopProp2DAcoVTIDenQ_DEO2_FDTD_stats(kwargs[:stats], kwargs[:ginsu], ntmod, time()-time1, cumtime_io, cumtime_ex, cumtime_im)
+    JopProp2DAcoVTIDenQ_DEO2_FDTD_stats(kwargs[:stats], kwargs[:ginsu], ntmod, time()-time1,
+        cumtime_io, cumtime_ex, cumtime_im)
 
     # undo ginsu
     for prop in keys(kwargs[:active_modelset])
@@ -904,7 +910,7 @@ end
 
 function Base.close(jet::Jet{D,R,typeof(JopProp2DAcoVTIDenQ_DEO2_FDTD_f!)}) where {D,R}
     rm("$(state(jet).srcfieldfile)-pold", force=true)
-    rm("$(state(jet).srcfieldfile)-mold", force=true)
     rm("$(state(jet).srcfieldfile)-pspace", force=true)
+    rm("$(state(jet).srcfieldfile)-mold", force=true)
     rm("$(state(jet).srcfieldfile)-mspace", force=true)
-end
\ No newline at end of file
+end
diff --git a/src/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl b/src/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl
index da06b3f..3e24a22 100644
--- a/src/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl
+++ b/src/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl
@@ -220,6 +220,9 @@ model assigned to it.
 The model is 4D with the slowest dimension for model parameter.  The earth model property maps to the slow-
 fimension array index via the `state(F).active_modelset` dictionary.
 
+The model is 4D with the slowest dimension for model parameter velocity [1] or buoyancy [2] in the case that both 
+`v` and `b` are active parameters. 
+        
 # Examples
 
 ## Model and acquisition geometry setup
@@ -445,7 +448,7 @@ function JopProp3DAcoIsoDenQ_DEO2_FDTD_nonlinearforward!(d::AbstractArray, m::Ab
 
     it0, ntmod_wav = WaveFD.default_ntmod(kwargs[:dtrec], kwargs[:dtmod], kwargs[:st], kwargs[:ntrec])
 
-    # source wavelet for injection, one for each source location
+    # source wavelet for injection, one for each source location 
     wavelet_realization = realizewavelet(kwargs[:wavelet], kwargs[:sz], kwargs[:sx], kwargs[:st], kwargs[:dtmod], ntmod_wav)
 
     # Get source and receiver interpolation coefficients
diff --git a/src/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl b/src/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl
index 80fd3d4..74c8794 100644
--- a/src/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl
+++ b/src/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl
@@ -1,10 +1,10 @@
 function JetProp3DAcoTTIDenQ_DEO2_FDTD(;
-        v = Float32[],
-        b = Float32[],
-        ϵ = Float32[],
-        η = Float32[],
-        θ = Float32[],
-        ϕ = Float32[],
+        v::Array{Float32,3} = ones(Float32,0,0,0),
+        b::Array{Float32,3} = ones(Float32,0,0,0),
+        ϵ::Array{Float32,3} = ones(Float32,0,0,0),
+        η::Array{Float32,3} = ones(Float32,0,0,0),
+        θ::Array{Float32,3} = ones(Float32,0,0,0),
+        ϕ::Array{Float32,3} = ones(Float32,0,0,0),
         f = 0.85,
         srcfieldfile = joinpath(tempdir(), "field-$(uuid4()).bin"),
         comptype = nothing,
@@ -57,7 +57,7 @@ function JetProp3DAcoTTIDenQ_DEO2_FDTD(;
         if length(x) > 0
             passive_modelset[n] = x
         else
-            active_modelset[n] = i # so that m[:,:,i] corresponds to property n∈("v","b","ϵ","η")
+            active_modelset[n] = i # so that m[:,:,:,i] corresponds to property n∈("v","b","ϵ","η")
             i += 1
         end
     end
@@ -119,7 +119,7 @@ function JetProp3DAcoTTIDenQ_DEO2_FDTD(;
     C = WaveFD.comptype(comptype, Float32)[1]
     compressor = Dict{String,WaveFD.Compressor{Float32,Float32,C,3}}()
 
-    # we need to serialize for the data, but not for the linearization
+    # we need to serialize "pold" for the data but, not (necessarily) for the linearization
     _active_wavefields = "pold" ∈ active_wavefields ? active_wavefields : [active_wavefields;"pold"]
 
     for active_wavefield in _active_wavefields
@@ -137,7 +137,6 @@ function JetProp3DAcoTTIDenQ_DEO2_FDTD(;
         error("Supplied imaging condition 'imgcondition' is not in [standard, FWI, RTM]")
     end
         
-    # construct:
     Jet(
         dom = dom,
         rng = rng,
@@ -371,7 +370,7 @@ Defaults for arguments are shown inside square brackets.
 * `θ [zeros(Float32,0,0,0)]` the symmetry axies tilt angle in radians. Assumed to be zero if not 
     specified. θ must be an array the same size as buoyancy `b`. 
 * `f [0.85]` See the discussion above regarding the **Pseudo-acoustic approximation** used here
-* `srcfieldfile ["field-$(uuid4()).bin"]` the full path to a scratch file used for 
+* `srcfieldfile [joinpath(tempdir(), "field-$(uuid4()).bin")]` the full path to a scratch file used for 
     the serializationof the compressed nonlinear source wavefield. 
 * `comptype [nothing]` the type of compression to use for the serialization of
     the nonlinear source wavefield. The type of compression must be one of:
@@ -412,16 +411,16 @@ Defaults for arguments are shown inside square brackets.
 * `dy [10.0]` Spacing of physical coordinates in the Y dimension.
 * `dx [10.0]` Spacing of physical coordinates in the X dimension.
 * `freqQ [5.0]` The center frequency for the Maxwell body approximation to dissipation only attenuation.
-    Please see `JetPackWaveFDFDFD` package documentation for more information concerning the attenuation model.   
+    Please see `JetPackWaveFD` package documentation for more information concerning the attenuation model.   
 * `qMin [0.1]` The minimum value for Qp at the boundary of the model used in our Maxwell body approximation
     to dissipation only attenuation. This is not a physically meaningful value for Qp, as we use the 
     attenuation to implement absorbing boundary conditions and eliminate outgoing waves on the 
     boundaries of the computational domain. 
-    Please see `JetPackWaveFDFDFD` package documentation for more information concerning the attenuation model.   
+    Please see `JetPackWaveFD` package documentation for more information concerning the attenuation model.   
 * `qInterior [100.0]` the value for Qp in the interior of the model used in our Maxwell body approximation
     to dissipation only attenuation. This is the value for Qp away from the absorbing boundaries and is 
     a physically meaningful value.
-    Please see `JetPackWaveFDFDFD` package documentation for more information concerning the attenuation model.   
+    Please see `JetPackWaveFD` package documentation for more information concerning the attenuation model.   
 * `padz, pady, padx [0.0], [0.0], [0.0]` - apply extra padding to the survey determined aperture in `Ginsu`.
     Please see `Ginsu` for more information. 
 * `nbz_cache, nby_cache, nbx_cache [512], [8], [8]` The size of cache blocks in the Z, X, and Y dimensions. 
@@ -477,10 +476,10 @@ function JopProp3DAcoTTIDenQ_DEO2_FDTD_nonlinearforward!(d::AbstractArray, m::Ab
         qMin = kwargs[:qMin],
         qInterior = kwargs[:qInterior])
 
-    # wave-fields
+    # wavefields
     wavefields = Dict("pcur"=>WaveFD.PCur(p), "pold"=>WaveFD.POld(p), "pspace"=>WaveFD.PSpace(p), "mcur"=>WaveFD.MCur(p), "mold"=>WaveFD.MOld(p), "mspace"=>WaveFD.MSpace(p))
 
-    # earth model
+    # ginsu'd earth model
     model_ginsu = Dict("v"=>WaveFD.V(p), "b"=>WaveFD.B(p), "ϵ"=>WaveFD.Eps(p), "η"=>WaveFD.Eta(p), "f"=>WaveFD.F(p), 
         "cosθ"=>WaveFD.CosTheta(p), "sinθ"=>WaveFD.SinTheta(p), "cosϕ"=>WaveFD.CosPhi(p), "sinϕ"=>WaveFD.SinPhi(p))
 
@@ -488,9 +487,12 @@ function JopProp3DAcoTTIDenQ_DEO2_FDTD_nonlinearforward!(d::AbstractArray, m::Ab
     for prop in keys(kwargs[:active_modelset])
         sub!(model_ginsu[prop], kwargs[:ginsu], @view(m[:,:,:,kwargs[:active_modelset][prop]]), extend=true)
     end
+    # ginsu'd earth model (passive-set)
     for prop in keys(kwargs[:passive_modelset])
         sub!(model_ginsu[prop], kwargs[:ginsu], kwargs[:passive_modelset][prop], extend=true)
     end
+
+    # TODO why no ginsu for f?
     model_ginsu["f"] .= kwargs[:f]
 
     ginsu_interior_range = interior(kwargs[:ginsu])
@@ -538,10 +540,13 @@ function JopProp3DAcoTTIDenQ_DEO2_FDTD_nonlinearforward!(d::AbstractArray, m::Ab
     if kwargs[:srcfieldfile] != ""
         for active_wavefield in active_wavefields
             filename = "$(kwargs[:srcfieldfile])-$(active_wavefield)"
-            if isfile(filename) == true
-                rm(filename)
+            try
+                isfile(filename) && rm(filename)
+                iofield[active_wavefield] = open(filename, "w")
+            catch
+                @info "Unable to open $(filename) on proc $(myid()) - $(gethostname())"
+                rethrow()
             end
-            iofield[active_wavefield] = open(filename, "w")
             open(kwargs[:compressor][active_wavefield])
         end
     end
@@ -554,7 +559,8 @@ function JopProp3DAcoTTIDenQ_DEO2_FDTD_nonlinearforward!(d::AbstractArray, m::Ab
     set_zero_subnormals(true)
     for it = 1:ntmod_wav
         if kwargs[:reportinterval] != 0 && (it % kwargs[:reportinterval] == 0 || it == ntmod_wav)
-            JopProp3DAcoTTIDenQ_DEO2_FDTD_write_history_nl(kwargs[:ginsu], it, ntmod_wav, time()-time1, cumtime_io, cumtime_ex, cumtime_pr, wavefields["pcur"], d)
+            JopProp3DAcoTTIDenQ_DEO2_FDTD_write_history_nl(kwargs[:ginsu], it, ntmod_wav, time()-time1,
+                cumtime_io, cumtime_ex, cumtime_pr, wavefields["pcur"], d)
         end
 
         # propagate and wavefield swap
@@ -580,11 +586,11 @@ function JopProp3DAcoTTIDenQ_DEO2_FDTD_nonlinearforward!(d::AbstractArray, m::Ab
             if kwargs[:srcfieldfile] != ""
                 cumtime_io += @elapsed for active_wavefield in active_wavefields
                     if kwargs[:isinterior]
-                        WaveFD.compressedwrite(iofield[active_wavefield], kwargs[:compressor][active_wavefield], div(it-1,itskip)+1,
-                            wavefields[active_wavefield], ginsu_interior_range)
+                        WaveFD.compressedwrite(iofield[active_wavefield], kwargs[:compressor][active_wavefield],
+                            div(it-1,itskip)+1, wavefields[active_wavefield], ginsu_interior_range)
                     else
-                        WaveFD.compressedwrite(iofield[active_wavefield], kwargs[:compressor][active_wavefield], div(it-1,itskip)+1,
-                            wavefields[active_wavefield])
+                        WaveFD.compressedwrite(iofield[active_wavefield], kwargs[:compressor][active_wavefield],
+                            div(it-1,itskip)+1, wavefields[active_wavefield])
                     end
                 end
             end
@@ -658,21 +664,22 @@ function JopProp3DAcoTTIDenQ_DEO2_FDTD_df!(δd::AbstractArray, δm::AbstractArra
         qMin = kwargs[:qMin],
         qInterior = kwargs[:qInterior])
 
-    # wave-fields
+    # wavefields
     pcur,pold = WaveFD.PCur(p),WaveFD.POld(p)
 
-    # earth model
-    earth_ginsu = Dict("v"=>WaveFD.V(p), "b"=>WaveFD.B(p), "ϵ"=>WaveFD.Eps(p), "η"=>WaveFD.Eta(p), "f"=>WaveFD.F(p), 
+    # ginsu'd earth model
+    model_ginsu = Dict("v"=>WaveFD.V(p), "b"=>WaveFD.B(p), "ϵ"=>WaveFD.Eps(p), "η"=>WaveFD.Eta(p), "f"=>WaveFD.F(p), 
         "cosθ"=>WaveFD.CosTheta(p), "sinθ"=>WaveFD.SinTheta(p), "cosϕ"=>WaveFD.CosPhi(p), "sinϕ"=>WaveFD.SinPhi(p))
 
     # ginsu'd earth model (active-set)
     for prop in keys(kwargs[:active_modelset])
-        sub!(earth_ginsu[prop], kwargs[:ginsu], @view(kwargs[:mₒ][:,:,:,kwargs[:active_modelset][prop]]), extend=true)
+        sub!(model_ginsu[prop], kwargs[:ginsu], @view(kwargs[:mₒ][:,:,:,kwargs[:active_modelset][prop]]), extend=true)
     end
+    # ginsu'd earth model (passive-set)
     for prop in keys(kwargs[:passive_modelset])
-        sub!(earth_ginsu[prop], kwargs[:ginsu], kwargs[:passive_modelset][prop], extend=true)
+        sub!(model_ginsu[prop], kwargs[:ginsu], kwargs[:passive_modelset][prop], extend=true)
     end
-    earth_ginsu["f"] .= kwargs[:f]
+    model_ginsu["f"] .= kwargs[:f]
 
     δm_ginsu = Dict{String,Array{Float32,3}}()
     for prop in keys(kwargs[:active_modelset])
@@ -682,7 +689,7 @@ function JopProp3DAcoTTIDenQ_DEO2_FDTD_df!(δd::AbstractArray, δm::AbstractArra
     ginsu_interior_range = interior(kwargs[:ginsu])
 
     # pre-compute receiver interpolation coefficients
-    local iz,iy,ix,c
+    local iz, iy, ix, c
     if kwargs[:interpmethod] == :hicks
         iz, iy, ix, c = WaveFD.hickscoeffs(kwargs[:dz], kwargs[:dy], kwargs[:dx], z0_ginsu, y0_ginsu, x0_ginsu, nz_ginsu, ny_ginsu, nx_ginsu, kwargs[:rz], kwargs[:ry], kwargs[:rx])
     else
@@ -730,11 +737,11 @@ function JopProp3DAcoTTIDenQ_DEO2_FDTD_df!(δd::AbstractArray, δm::AbstractArra
             # read source field from disk
             cumtime_io += @elapsed for active_wavefield in kwargs[:active_wavefields]
                 if kwargs[:isinterior]
-                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield], div(it-1,itskip)+1,
-                        wavefields[active_wavefield], ginsu_interior_range)
+                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield],
+                        div(it-1,itskip)+1, wavefields[active_wavefield], ginsu_interior_range)
                 else
-                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield], div(it-1,itskip)+1,
-                        wavefields[active_wavefield])
+                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield],
+                        div(it-1,itskip)+1, wavefields[active_wavefield])
                 end
             end
 
@@ -783,22 +790,22 @@ function JopProp3DAcoTTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
         qMin = kwargs[:qMin],
         qInterior = kwargs[:qInterior])
 
-    # wave-fields
+    # wavefields
     pcur,pold = WaveFD.PCur(p),WaveFD.POld(p)
 
     # ginsu'd earth model
-    earth_ginsu = Dict("v"=>WaveFD.V(p), "b"=>WaveFD.B(p), "ϵ"=>WaveFD.Eps(p), "η"=>WaveFD.Eta(p), "f"=>WaveFD.F(p), 
+    model_ginsu = Dict("v"=>WaveFD.V(p), "b"=>WaveFD.B(p), "ϵ"=>WaveFD.Eps(p), "η"=>WaveFD.Eta(p), "f"=>WaveFD.F(p), 
         "cosθ"=>WaveFD.CosTheta(p), "sinθ"=>WaveFD.SinTheta(p), "cosϕ"=>WaveFD.CosPhi(p), "sinϕ"=>WaveFD.SinPhi(p))
 
-    # active model-set
+    # ginsu'd earth model (active-set)
     for prop in keys(kwargs[:active_modelset])
-        sub!(earth_ginsu[prop], kwargs[:ginsu], @view(kwargs[:mₒ][:,:,:,kwargs[:active_modelset][prop]]), extend=true)
+        sub!(model_ginsu[prop], kwargs[:ginsu], @view(kwargs[:mₒ][:,:,:,kwargs[:active_modelset][prop]]), extend=true)
     end
-    # passive model-set
+    # ginsu'd earth model (passive-set)
     for prop in keys(kwargs[:passive_modelset])
-        sub!(earth_ginsu[prop], kwargs[:ginsu], kwargs[:passive_modelset][prop], extend=true)
+        sub!(model_ginsu[prop], kwargs[:ginsu], kwargs[:passive_modelset][prop], extend=true)
     end
-    earth_ginsu["f"] .= kwargs[:f]
+    model_ginsu["f"] .= kwargs[:f]
 
     δm_ginsu = Dict{String,Array{Float32,3}}()
     for prop in keys(kwargs[:active_modelset])
@@ -808,7 +815,7 @@ function JopProp3DAcoTTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
     ginsu_interior_range = interior(kwargs[:ginsu])
 
     # Get receiver interpolation coefficients
-    local iz,iy,ix,c
+    local iz, iy, ix, c
     if kwargs[:interpmethod] == :hicks
         iz, iy, ix, c = WaveFD.hickscoeffs(kwargs[:dz], kwargs[:dy], kwargs[:dx], z0_ginsu, y0_ginsu, x0_ginsu, 
             nz_ginsu, ny_ginsu, nx_ginsu, kwargs[:rz], kwargs[:ry], kwargs[:rx])
@@ -821,7 +828,7 @@ function JopProp3DAcoTTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
             kz = iz[i][jz]
             ky = iy[i][jy]
             kx = ix[i][jx]
-            c[i][jz,jy,jx] *= kwargs[:dtmod]^2 * earth_ginsu["v"][kz,ky,kx]^2 / earth_ginsu["b"][kz,ky,kx]
+            c[i][jz,jy,jx] *= kwargs[:dtmod]^2 * model_ginsu["v"][kz,ky,kx]^2 / model_ginsu["b"][kz,ky,kx]
         end
     end
     blks = WaveFD.source_blocking(nz_ginsu, ny_ginsu, nx_ginsu, kwargs[:nbz_inject], kwargs[:nby_inject], kwargs[:nbx_inject], iz, iy, ix, c)
@@ -857,7 +864,8 @@ function JopProp3DAcoTTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
     set_zero_subnormals(true)
     for it = ntmod:-1:1
         if kwargs[:reportinterval] != 0 && (it % kwargs[:reportinterval] == 0 || it == ntmod)
-            JopProp3DAcoTTIDenQ_DEO2_FDTD_write_history_ln(kwargs[:ginsu], it, ntmod, time()-time1, cumtime_io, cumtime_ex, cumtime_im, cumtime_pr, pcur, δdinterp, "adjoint")
+            JopProp3DAcoTTIDenQ_DEO2_FDTD_write_history_ln(kwargs[:ginsu], it, ntmod, time()-time1,
+                cumtime_io, cumtime_ex, cumtime_im, cumtime_pr, pcur, δdinterp, "adjoint")
         end
 
         # propagate and wavefield swap
@@ -871,11 +879,11 @@ function JopProp3DAcoTTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
             # read source field from disk
             cumtime_io += @elapsed for active_wavefield in kwargs[:active_wavefields]
                 if kwargs[:isinterior]
-                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield], div(it-1,itskip)+1,
-                        wavefields[active_wavefield], ginsu_interior_range)
+                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield],
+                        div(it-1,itskip)+1, wavefields[active_wavefield], ginsu_interior_range)
                 else
-                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield], div(it-1,itskip)+1,
-                        wavefields[active_wavefield])
+                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield],
+                        div(it-1,itskip)+1, wavefields[active_wavefield])
                 end
             end
 
@@ -884,7 +892,8 @@ function JopProp3DAcoTTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
         end
     end
     set_zero_subnormals(false)
-    JopProp3DAcoTTIDenQ_DEO2_FDTD_stats(kwargs[:stats], kwargs[:ginsu], ntmod, time()-time1, cumtime_io, cumtime_ex, cumtime_im)
+    JopProp3DAcoTTIDenQ_DEO2_FDTD_stats(kwargs[:stats], kwargs[:ginsu], ntmod, time()-time1,
+        cumtime_io, cumtime_ex, cumtime_im)
 
     # undo ginsu
     for prop in keys(kwargs[:active_modelset])
@@ -901,11 +910,11 @@ function JopProp3DAcoTTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
     δm
 end
 
-modelindex(F::Jop{T}, key::AbstractString) where {D,R,T<:Jet{D,R,typeof(JopProp3DAcoTTIDenQ_DEO2_FDTD_f!)}} = state(F).active_modelset[key]
+modelindex(F::Jop{T}, key) where {D,R,T<:Jet{D,R,typeof(JopProp3DAcoTTIDenQ_DEO2_FDTD_f!)}} = state(F).active_modelset[key]
 
-function srcillum!(γ, A::T, m::AbstractArray{Float32}) where {D,R,J<:Jet{D,R,typeof(JopProp3DAcoTTIDenQ_DEO2_FDTD_f!)},T<:Jop{J}}
+function srcillum!(γ, A::Jop{T}, m::AbstractArray{Float32}) where {D,R,T<:Jet{D,R,typeof(JopProp3DAcoTTIDenQ_DEO2_FDTD_f!)}}
     s = state(A)
-    isvalid, _chksum = isvalid_srcfieldfile(m, s.srcfieldhost[], s.srcfieldfile*"-pold", s.chksum[])
+    isvalid, _chksum = isvalid_srcfieldfile(jet(A).mₒ, s.srcfieldhost[], s.srcfieldfile*"-pold", s.chksum[])
     if !isvalid
         JopProp3DAcoTTIDenQ_DEO2_FDTD_nonlinearforward!(Array{Float32}(undef,0,0), m; s...)
         s.chksum[] = _chksum
@@ -954,9 +963,9 @@ end
     @info "Prop3DAcoTTIDenQ_DEO2_FDTD, nonlinear forward, time step $kt of $nt $mcells MCells/s (IO=$IO, EX=$EX, PR=$PR) -- rms d,p; $rmsd $rmsp"
 end
 
-function Base.close(j::Jet{D,R,typeof(JopProp3DAcoTTIDenQ_DEO2_FDTD_f!)}) where {D,R}
-    rm("$(state(j).srcfieldfile)-pold", force=true)
-    rm("$(state(j).srcfieldfile)-mold", force=true)
-    rm("$(state(j).srcfieldfile)-pspace", force=true)
-    rm("$(state(j).srcfieldfile)-mspace", force=true)
-end
\ No newline at end of file
+function Base.close(jet::Jet{D,R,typeof(JopProp3DAcoTTIDenQ_DEO2_FDTD_f!)}) where {D,R}
+    rm("$(state(jet).srcfieldfile)-pold", force=true)
+    rm("$(state(jet).srcfieldfile)-pspace", force=true)
+    rm("$(state(jet).srcfieldfile)-mold", force=true)
+    rm("$(state(jet).srcfieldfile)-mspace", force=true)
+end
diff --git a/src/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl b/src/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl
index 8e46bc3..0436bd0 100644
--- a/src/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl
+++ b/src/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl
@@ -1,8 +1,8 @@
 function JetProp3DAcoVTIDenQ_DEO2_FDTD(;
-        v = Float32[],
-        b = Float32[],
-        ϵ = Float32[],
-        η = Float32[],
+        v::Array{Float32,3} = ones(Float32,0,0,0),
+        b::Array{Float32,3} = ones(Float32,0,0,0),
+        ϵ::Array{Float32,3} = ones(Float32,0,0,0),
+        η::Array{Float32,3} = ones(Float32,0,0,0),
         f = 0.85,
         srcfieldfile = joinpath(tempdir(), "field-$(uuid4()).bin"),
         comptype = nothing,
@@ -55,12 +55,14 @@ function JetProp3DAcoVTIDenQ_DEO2_FDTD(;
         if length(x) > 0
             passive_modelset[n] = x
         else
-            active_modelset[n] = i # so that m[:,:,i] corresponds to property n∈("v","b","ϵ","η")
+            active_modelset[n] = i # so that m[:,:,:,i] corresponds to property n∈("v","b","ϵ","η")
             i += 1
         end
     end
     @assert length(active_modelset) > 0
 
+    # TODO do we need some assert for sensible buoyancy?
+
     # active and passive wavefields (an active wavefield is serialized to disk)
     active_modelset_keys = keys(active_modelset)
     local active_wavefields, modeltype
@@ -113,7 +115,7 @@ function JetProp3DAcoVTIDenQ_DEO2_FDTD(;
     C = WaveFD.comptype(comptype, Float32)[1]
     compressor = Dict{String,WaveFD.Compressor{Float32,Float32,C,3}}()
 
-    # we need to serialize for the data, but not for the linearization
+    # we need to serialize "pold" for the data but, not (necessarily) for the linearization
     _active_wavefields = "pold" ∈ active_wavefields ? active_wavefields : [active_wavefields;"pold"]
 
     for active_wavefield in _active_wavefields
@@ -131,7 +133,6 @@ function JetProp3DAcoVTIDenQ_DEO2_FDTD(;
         error("Supplied imaging condition 'imgcondition' is not in [standard, FWI, RTM]")
     end
 
-    # construct:
     Jet(
         dom = dom,
         rng = rng,
@@ -339,7 +340,7 @@ J = JopLnProp3DAcoVTIDenQ_DEO2_FDTD(; m₀ = m₀, b = b, isinterior=true, nspon
 
 # Required Parameters
 * `m₀` the point at which the jet is linearized. Note this argument is required in the constuctor for  
-    `JopLnProp2DAcoVTIDenQ_DEO2_FDTD` but not `JopNlProp2DAcoVTIDenQ_DEO2_FDTD`. This constuctor is shown  
+    `JopLnProp3DAcoVTIDenQ_DEO2_FDTD` but not `JopNlProp3DAcoVTIDenQ_DEO2_FDTD`. This constuctor is shown  
     in the last example above. Please note that you must consider which parameters are active and passive,
     per the discussion on model parameters above and examples.
 * `dtmod` The sample rate for the modeled data. You can establish a lower bound for the modeling sample rate  
@@ -463,19 +464,22 @@ function JopProp3DAcoVTIDenQ_DEO2_FDTD_nonlinearforward!(d::AbstractArray, m::Ab
         qMin = kwargs[:qMin],
         qInterior = kwargs[:qInterior])
 
-    # wave-fields
+    # wavefields
     wavefields = Dict("pcur"=>WaveFD.PCur(p), "pold"=>WaveFD.POld(p), "pspace"=>WaveFD.PSpace(p), "mcur"=>WaveFD.MCur(p), "mold"=>WaveFD.MOld(p), "mspace"=>WaveFD.MSpace(p))
 
-    # earth model
+    # ginsu'd earth model
     model_ginsu = Dict("v"=>WaveFD.V(p), "b"=>WaveFD.B(p), "ϵ"=>WaveFD.Eps(p), "η"=>WaveFD.Eta(p), "f"=>WaveFD.F(p))
 
     # ginsu'd earth model (active-set)
     for prop in keys(kwargs[:active_modelset])
         sub!(model_ginsu[prop], kwargs[:ginsu], @view(m[:,:,:,kwargs[:active_modelset][prop]]), extend=true)
     end
+    # ginsu'd earth model (passive-set)
     for prop in keys(kwargs[:passive_modelset])
         sub!(model_ginsu[prop], kwargs[:ginsu], kwargs[:passive_modelset][prop], extend=true)
     end
+
+    # TODO why no ginsu for f?
     model_ginsu["f"] .= kwargs[:f]
 
     ginsu_interior_range = interior(kwargs[:ginsu])
@@ -523,10 +527,13 @@ function JopProp3DAcoVTIDenQ_DEO2_FDTD_nonlinearforward!(d::AbstractArray, m::Ab
     if kwargs[:srcfieldfile] != ""
         for active_wavefield in active_wavefields
             filename = "$(kwargs[:srcfieldfile])-$(active_wavefield)"
-            if isfile(filename) == true
-                rm(filename)
+            try
+                isfile(filename) && rm(filename)
+                iofield[active_wavefield] = open(filename, "w")
+            catch
+                @info "Unable to open $(filename) on proc $(myid()) - $(gethostname())"
+                rethrow()
             end
-            iofield[active_wavefield] = open(filename, "w")
             open(kwargs[:compressor][active_wavefield])
         end
     end
@@ -539,7 +546,8 @@ function JopProp3DAcoVTIDenQ_DEO2_FDTD_nonlinearforward!(d::AbstractArray, m::Ab
     set_zero_subnormals(true)
     for it = 1:ntmod_wav
         if kwargs[:reportinterval] != 0 && (it % kwargs[:reportinterval] == 0 || it == ntmod_wav)
-            JopProp3DAcoVTIDenQ_DEO2_FDTD_write_history_nl(kwargs[:ginsu], it, ntmod_wav, time()-time1, cumtime_io, cumtime_ex, cumtime_pr, wavefields["pcur"], d)
+            JopProp3DAcoVTIDenQ_DEO2_FDTD_write_history_nl(kwargs[:ginsu], it, ntmod_wav, time()-time1,
+                cumtime_io, cumtime_ex, cumtime_pr, wavefields["pcur"], d)
         end
 
         # propagate and wavefield swap
@@ -565,11 +573,11 @@ function JopProp3DAcoVTIDenQ_DEO2_FDTD_nonlinearforward!(d::AbstractArray, m::Ab
             if kwargs[:srcfieldfile] != ""
                 cumtime_io += @elapsed for active_wavefield in active_wavefields
                     if kwargs[:isinterior]
-                        WaveFD.compressedwrite(iofield[active_wavefield], kwargs[:compressor][active_wavefield], div(it-1,itskip)+1,
-                            wavefields[active_wavefield], ginsu_interior_range)
+                        WaveFD.compressedwrite(iofield[active_wavefield], kwargs[:compressor][active_wavefield],
+                            div(it-1,itskip)+1, wavefields[active_wavefield], ginsu_interior_range)
                     else
-                        WaveFD.compressedwrite(iofield[active_wavefield], kwargs[:compressor][active_wavefield], div(it-1,itskip)+1,
-                            wavefields[active_wavefield])
+                        WaveFD.compressedwrite(iofield[active_wavefield], kwargs[:compressor][active_wavefield],
+                            div(it-1,itskip)+1, wavefields[active_wavefield])
                     end
                 end
             end
@@ -643,20 +651,21 @@ function JopProp3DAcoVTIDenQ_DEO2_FDTD_df!(δd::AbstractArray, δm::AbstractArra
         qMin = kwargs[:qMin],
         qInterior = kwargs[:qInterior])
 
-    # wave-fields
+    # wavefields
     pcur,pold = WaveFD.PCur(p),WaveFD.POld(p)
 
-    # earth model
-    earth_ginsu = Dict("v"=>WaveFD.V(p), "b"=>WaveFD.B(p), "ϵ"=>WaveFD.Eps(p), "η"=>WaveFD.Eta(p), "f"=>WaveFD.F(p))
+    # ginsu'd earth model
+    model_ginsu = Dict("v"=>WaveFD.V(p), "b"=>WaveFD.B(p), "ϵ"=>WaveFD.Eps(p), "η"=>WaveFD.Eta(p), "f"=>WaveFD.F(p))
 
     # ginsu'd earth model (active-set)
     for prop in keys(kwargs[:active_modelset])
-        sub!(earth_ginsu[prop], kwargs[:ginsu], @view(kwargs[:mₒ][:,:,:,kwargs[:active_modelset][prop]]), extend=true)
+        sub!(model_ginsu[prop], kwargs[:ginsu], @view(kwargs[:mₒ][:,:,:,kwargs[:active_modelset][prop]]), extend=true)
     end
+    # ginsu'd earth model (passive-set)
     for prop in keys(kwargs[:passive_modelset])
-        sub!(earth_ginsu[prop], kwargs[:ginsu], kwargs[:passive_modelset][prop], extend=true)
+        sub!(model_ginsu[prop], kwargs[:ginsu], kwargs[:passive_modelset][prop], extend=true)
     end
-    earth_ginsu["f"] .= kwargs[:f]
+    model_ginsu["f"] .= kwargs[:f]
 
     δm_ginsu = Dict{String,Array{Float32,3}}()
     for prop in keys(kwargs[:active_modelset])
@@ -666,7 +675,7 @@ function JopProp3DAcoVTIDenQ_DEO2_FDTD_df!(δd::AbstractArray, δm::AbstractArra
     ginsu_interior_range = interior(kwargs[:ginsu])
 
     # pre-compute receiver interpolation coefficients
-    local iz,iy,ix,c
+    local iz, iy, ix, c
     if kwargs[:interpmethod] == :hicks
         iz, iy, ix, c = WaveFD.hickscoeffs(kwargs[:dz], kwargs[:dy], kwargs[:dx], z0_ginsu, y0_ginsu, x0_ginsu, nz_ginsu, ny_ginsu, nx_ginsu, kwargs[:rz], kwargs[:ry], kwargs[:rx])
     else
@@ -714,11 +723,11 @@ function JopProp3DAcoVTIDenQ_DEO2_FDTD_df!(δd::AbstractArray, δm::AbstractArra
             # read source field from disk
             cumtime_io += @elapsed for active_wavefield in kwargs[:active_wavefields]
                 if kwargs[:isinterior]
-                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield], div(it-1,itskip)+1,
-                        wavefields[active_wavefield], ginsu_interior_range)
+                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield],
+                        div(it-1,itskip)+1, wavefields[active_wavefield], ginsu_interior_range)
                 else
-                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield], div(it-1,itskip)+1,
-                        wavefields[active_wavefield])
+                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield],
+                        div(it-1,itskip)+1, wavefields[active_wavefield])
                 end
             end
 
@@ -767,21 +776,21 @@ function JopProp3DAcoVTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
         qMin = kwargs[:qMin],
         qInterior = kwargs[:qInterior])
 
-    # wave-fields
+    # wavefields
     pcur,pold = WaveFD.PCur(p),WaveFD.POld(p)
 
     # ginsu'd earth model
-    earth_ginsu = Dict("v"=>WaveFD.V(p), "b"=>WaveFD.B(p), "ϵ"=>WaveFD.Eps(p), "η"=>WaveFD.Eta(p), "f"=>WaveFD.F(p))
+    model_ginsu = Dict("v"=>WaveFD.V(p), "b"=>WaveFD.B(p), "ϵ"=>WaveFD.Eps(p), "η"=>WaveFD.Eta(p), "f"=>WaveFD.F(p))
 
-    # active model-set
+    # ginsu'd earth model (active-set)
     for prop in keys(kwargs[:active_modelset])
-        sub!(earth_ginsu[prop], kwargs[:ginsu], @view(kwargs[:mₒ][:,:,:,kwargs[:active_modelset][prop]]), extend=true)
+        sub!(model_ginsu[prop], kwargs[:ginsu], @view(kwargs[:mₒ][:,:,:,kwargs[:active_modelset][prop]]), extend=true)
     end
-    # passive model-set
+    # ginsu'd earth model (passive-set)
     for prop in keys(kwargs[:passive_modelset])
-        sub!(earth_ginsu[prop], kwargs[:ginsu], kwargs[:passive_modelset][prop], extend=true)
+        sub!(model_ginsu[prop], kwargs[:ginsu], kwargs[:passive_modelset][prop], extend=true)
     end
-    earth_ginsu["f"] .= kwargs[:f]
+    model_ginsu["f"] .= kwargs[:f]
 
     δm_ginsu = Dict{String,Array{Float32,3}}()
     for prop in keys(kwargs[:active_modelset])
@@ -804,7 +813,7 @@ function JopProp3DAcoVTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
             kz = iz[i][jz]
             ky = iy[i][jy]
             kx = ix[i][jx]
-            c[i][jz,jy,jx] *= kwargs[:dtmod]^2 * earth_ginsu["v"][kz,ky,kx]^2 / earth_ginsu["b"][kz,ky,kx]
+            c[i][jz,jy,jx] *= kwargs[:dtmod]^2 * model_ginsu["v"][kz,ky,kx]^2 / model_ginsu["b"][kz,ky,kx]
         end
     end
     blks = WaveFD.source_blocking(nz_ginsu, ny_ginsu, nx_ginsu, kwargs[:nbz_inject], kwargs[:nby_inject], kwargs[:nbx_inject], iz, iy, ix, c)
@@ -840,7 +849,8 @@ function JopProp3DAcoVTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
     set_zero_subnormals(true)
     for it = ntmod:-1:1
         if kwargs[:reportinterval] != 0 && (it % kwargs[:reportinterval] == 0 || it == ntmod)
-            JopProp3DAcoVTIDenQ_DEO2_FDTD_write_history_ln(kwargs[:ginsu], it, ntmod, time()-time1, cumtime_io, cumtime_ex, cumtime_im, cumtime_pr, pcur, δdinterp, "adjoint")
+            JopProp3DAcoVTIDenQ_DEO2_FDTD_write_history_ln(kwargs[:ginsu], it, ntmod, time()-time1,
+                cumtime_io, cumtime_ex, cumtime_im, cumtime_pr, pcur, δdinterp, "adjoint")
         end
 
         # propagate and wavefield swap
@@ -854,11 +864,11 @@ function JopProp3DAcoVTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
             # read source field from disk
             cumtime_io += @elapsed for active_wavefield in kwargs[:active_wavefields]
                 if kwargs[:isinterior]
-                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield], div(it-1,itskip)+1,
-                        wavefields[active_wavefield], ginsu_interior_range)
+                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield],
+                        div(it-1,itskip)+1, wavefields[active_wavefield], ginsu_interior_range)
                 else
-                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield], div(it-1,itskip)+1,
-                        wavefields[active_wavefield])
+                    WaveFD.compressedread!(iofields[active_wavefield], kwargs[:compressor][active_wavefield],
+                        div(it-1,itskip)+1, wavefields[active_wavefield])
                 end
             end
 
@@ -867,7 +877,8 @@ function JopProp3DAcoVTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
         end
     end
     set_zero_subnormals(false)
-    JopProp3DAcoVTIDenQ_DEO2_FDTD_stats(kwargs[:stats], kwargs[:ginsu], ntmod, time()-time1, cumtime_io, cumtime_ex, cumtime_im)
+    JopProp3DAcoVTIDenQ_DEO2_FDTD_stats(kwargs[:stats], kwargs[:ginsu], ntmod, time()-time1,
+        cumtime_io, cumtime_ex, cumtime_im)
 
     # undo ginsu
     for prop in keys(kwargs[:active_modelset])
@@ -884,11 +895,11 @@ function JopProp3DAcoVTIDenQ_DEO2_FDTD_df′!(δm::AbstractArray, δd::AbstractA
     δm
 end
 
-modelindex(F::Jop{T}, key::AbstractString) where {D,R,T<:Jet{D,R,typeof(JopProp3DAcoVTIDenQ_DEO2_FDTD_f!)}} = state(F).active_modelset[key]
+modelindex(F::Jop{T}, key) where {D,R,T<:Jet{D,R,typeof(JopProp3DAcoVTIDenQ_DEO2_FDTD_f!)}} = state(F).active_modelset[key]
 
-function srcillum!(γ, A::T, m::AbstractArray{Float32}) where {D,R,J<:Jet{D,R,typeof(JopProp3DAcoVTIDenQ_DEO2_FDTD_f!)},T<:Jop{J}}
+function srcillum!(γ, A::Jop{T}, m::AbstractArray{Float32}) where {D,R,T<:Jet{D,R,typeof(JopProp3DAcoVTIDenQ_DEO2_FDTD_f!)}}
     s = state(A)
-    isvalid, _chksum = isvalid_srcfieldfile(m, s.srcfieldhost[], s.srcfieldfile*"-pold", s.chksum[])
+    isvalid, _chksum = isvalid_srcfieldfile(jet(A).mₒ, s.srcfieldhost[], s.srcfieldfile*"-pold", s.chksum[])
     if !isvalid
         JopProp3DAcoVTIDenQ_DEO2_FDTD_nonlinearforward!(Array{Float32}(undef,0,0), m; s...)
         s.chksum[] = _chksum
@@ -937,9 +948,9 @@ end
     @info "Prop3DAcoVTIDenQ_DEO2_FDTD, nonlinear forward, time step $kt of $nt $mcells MCells/s (IO=$IO, EX=$EX, PR=$PR) -- rms d,p; $rmsd $rmsp"
 end
 
-function Base.close(j::Jet{D,R,typeof(JopProp3DAcoVTIDenQ_DEO2_FDTD_f!)}) where {D,R}
-    rm("$(state(j).srcfieldfile)-pold", force=true)
-    rm("$(state(j).srcfieldfile)-mold", force=true)
-    rm("$(state(j).srcfieldfile)-pspace", force=true)
-    rm("$(state(j).srcfieldfile)-mspace", force=true)
-end
\ No newline at end of file
+function Base.close(jet::Jet{D,R,typeof(JopProp3DAcoVTIDenQ_DEO2_FDTD_f!)}) where {D,R}
+    rm("$(state(jet).srcfieldfile)-pold", force=true)
+    rm("$(state(jet).srcfieldfile)-pspace", force=true)
+    rm("$(state(jet).srcfieldfile)-mold", force=true)
+    rm("$(state(jet).srcfieldfile)-mspace", force=true)
+end
diff --git a/test/jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl b/test/jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl
index 83ac569..d686540 100644
--- a/test/jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl
+++ b/test/jop_prop2DAcoIsoDenQ_DEO2_FDTD.jl
@@ -18,8 +18,22 @@ function make_op(interpmethod, modeltype, fs; comptype = Float32)
 
     wavelet = WaveletCausalRicker(f=5.0)
 
-    b = ones(Float32,nz,nx)
-    v = 1500 .* ones(Float32,nz,nx)
+    b = Float32(1)    .* ones(Float32,nz,nx)
+    v = Float32(1500) .* ones(Float32,nz,nx)
+
+    local m
+    if modeltype == WaveFD.Prop2DAcoIsoDenQ_DEO2_FDTD_Model_V
+        kwargs = (b = b, )
+        m = reshape(v, nz, nx, 1)
+    elseif modeltype == WaveFD.Prop2DAcoIsoDenQ_DEO2_FDTD_Model_VB
+        kwargs = (nz = nz, nx = nx)
+        m = zeros(Float32, nz, nx, 2)
+        m[:,:,1] .= v
+        m[:,:,2] .= b
+    elseif modeltype == WaveFD.Prop2DAcoIsoDenQ_DEO2_FDTD_Model_B
+        kwargs = (v = v, )
+        m = reshape(b, nz, nx, 1)
+    end
 
     local m
     if modeltype == WaveFD.Prop2DAcoIsoDenQ_DEO2_FDTD_Model_V
@@ -188,7 +202,7 @@ end
                 end
             end
             u_ana[:,iz,ix] = irfft(U,nt_pad)[1:ntrec]
-            printfmt("{:.2f}\r", ((ix-1)*nz+iz-1)/((nz-1)*(nx-1))*100))
+            printfmt("{:.2f}\r", ((ix-1)*nz+iz-1)/((nz-1)*(nx-1))*100)
         end
         u_ana
     end
diff --git a/test/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl b/test/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl
index 9c5b62d..e1c2dfa 100644
--- a/test/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl
+++ b/test/jop_prop2DAcoTTIDenQ_DEO2_FDTD.jl
@@ -153,8 +153,7 @@ end
     close(F)
 end
 
-#=
-@testset "JopProp2DAcoTTIDenQ_DEO2_FDTD -- analytic, direct, for model type $modeltype, interpmethod=$interpmethod" for modeltype in ("vϵη","v"), interpmethod in (:linear,:hicks)
+@test_skip @testset "JopProp2DAcoTTIDenQ_DEO2_FDTD -- analytic, direct, for model type $modeltype, interpmethod=$interpmethod" for modeltype in ("vϵη","v"), interpmethod in (:linear,:hicks)
     z,x,dz,dx,dtrec,dtmod,t,sz,sx,c = 0.5,1.0,0.02,0.02,0.004,0.0004,1.0,0.25,0.02,1.5
     nz,nx,ntrec = round(Int,z/dz)+1,round(Int,x/dx)+1,round(Int,t/dtrec)+1
 
@@ -184,7 +183,7 @@ end
                 end
             end
             u_ana[:,iz,ix] = irfft(U,nt_pad)[1:ntrec]
-            printfmt("{:.2f}\r", ((ix-1)*nz+iz-1)/((nz-1)*(nx-1))*100))
+            printfmt("{:.2f}\r", ((ix-1)*nz+iz-1)/((nz-1)*(nx-1))*100)
         end
         u_ana
     end
@@ -271,5 +270,4 @@ end
     @test cor(amp_fd[:],amp_ana[:]) > .99
 
     close(M)
-end
-=#
\ No newline at end of file
+end
\ No newline at end of file
diff --git a/test/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl b/test/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl
index 8088c73..2af9c55 100644
--- a/test/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl
+++ b/test/jop_prop2DAcoVTIDenQ_DEO2_FDTD.jl
@@ -156,8 +156,7 @@ end
     close(F)
 end
 
-#=
-@testset "JopProp2DAcoVTIDenQ_DEO2_FDTD -- analytic, direct, for model type $modeltype, interpmethod=$interpmethod" for modeltype in ("vϵη","v"), interpmethod in (:hicks,:linear)
+@test_skip @testset "JopProp2DAcoVTIDenQ_DEO2_FDTD -- analytic, direct, for model type $modeltype, interpmethod=$interpmethod" for modeltype in ("vϵη","v"), interpmethod in (:hicks,:linear)
     z,x,dz,dx,dtrec,dtmod,t,sz,sx,c = 0.5,1.0,0.02,0.02,0.004,0.0004,1.0,0.25,0.02,1.5
     nz,nx,ntrec = round(Int,z/dz)+1,round(Int,x/dx)+1,round(Int,t/dtrec)+1
 
@@ -187,7 +186,7 @@ end
                 end
             end
             u_ana[:,iz,ix] = irfft(U,nt_pad)[1:ntrec]
-            printfmt("{:.2f}\r", ((ix-1)*nz+iz-1)/((nz-1)*(nx-1))*100))
+            printfmt("{:.2f}\r", ((ix-1)*nz+iz-1)/((nz-1)*(nx-1))*100)
         end
         u_ana
     end
@@ -274,5 +273,4 @@ end
     @test cor(amp_fd[:],amp_ana[:]) > .99
 
     close(M)
-end
-=#
\ No newline at end of file
+end
\ No newline at end of file
diff --git a/test/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl b/test/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl
index cda94e3..3e75663 100644
--- a/test/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl
+++ b/test/jop_prop3DAcoIsoDenQ_DEO2_FDTD.jl
@@ -21,8 +21,22 @@ function make_op(interpmethod, modeltype, fs; comptype = Float32)
 
     wavelet = WaveletCausalRicker(f=5.0)
 
-    b = ones(Float32,nz,ny,nx)
-    v = 1500 .* ones(Float32,nz,ny,nx)
+    b = Float32(1)    .* ones(Float32,nz,ny,nx)
+    v = Float32(1500) .* ones(Float32,nz,ny,nx)
+
+    local m
+    if modeltype == WaveFD.Prop3DAcoIsoDenQ_DEO2_FDTD_Model_V
+        kwargs = (b = b, )
+        m = reshape(v, nz, ny, nx, 1)
+    elseif modeltype == WaveFD.Prop3DAcoIsoDenQ_DEO2_FDTD_Model_VB
+        kwargs = (nz = nz, ny = ny, nx = nx)
+        m = zeros(Float32, nz, ny, nx, 2)
+        m[:,:,:,1] .= v
+        m[:,:,:,2] .= b
+    elseif modeltype == WaveFD.Prop3DAcoIsoDenQ_DEO2_FDTD_Model_B
+        kwargs = (v = v, )
+        m = reshape(b, nz, ny, nx, 1)
+    end
 
     local m
     if modeltype == WaveFD.Prop3DAcoIsoDenQ_DEO2_FDTD_Model_V
@@ -161,8 +175,7 @@ end
     close(F)
 end
 
-#=
-@testset "JopProp3DAcoIsoDenQ_DEO2_FDTD -- analytic, direct, interpmethod=$interpmethod" for interpmethod in (:hicks,:linear)
+@test_skip @testset "JopProp3DAcoIsoDenQ_DEO2_FDTD -- analytic, direct, interpmethod=$interpmethod" for interpmethod in (:hicks,:linear)
     z,y,x,dz,dy,dx,dtrec,dtmod,tmax,sz,sy,sx,c = 2000.0,500.0,2500.0,40.0,40.0,40.0,0.016,0.004,2.0,1000.0,250.0,20.0,1500.0
     nz,ny,nx,nt = round(Int,z/dz)+1,round(Int,y/dy)+1,round(Int,x/dx)+1,round(Int,tmax/dtrec)+1
 
@@ -261,4 +274,5 @@ end
 
     close(M)
 end
-=#
+
+nothing
\ No newline at end of file
diff --git a/test/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl b/test/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl
index e6903de..484c940 100644
--- a/test/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl
+++ b/test/jop_prop3DAcoTTIDenQ_DEO2_FDTD.jl
@@ -157,8 +157,7 @@ end
     close(F)
 end
 
-#=
-@testset "JopProp3DAcoTTIDenQ_DEO2_FDTD -- analytic, direct, for model type $modeltype, interpmethod=$interpmethod" for modeltype in ("vϵη","v"), interpmethod in (:hicks,:linear)
+@test_skip @testset "JopProp3DAcoTTIDenQ_DEO2_FDTD -- analytic, direct, for model type $modeltype, interpmethod=$interpmethod" for modeltype in ("vϵη","v"), interpmethod in (:hicks,:linear)
     z,y,x,dz,dy,dx,dtrec,dtmod,tmax,sz,sy,sx,c = 2000.0,500.0,2500.0,40.0,40.0,40.0,0.016,0.004,2.0,1000.0,250.0,20.0,1500.0
     nz,ny,nx,nt = round(Int,z/dz)+1,round(Int,y/dy)+1,round(Int,x/dx)+1,round(Int,tmax/dtrec)+1
 
@@ -189,7 +188,7 @@ end
                 end
             end
             u_ana[:,iz,iy,ix] = irfft(U,nt_pad)[1:nt]
-            printfmt("{:.2f}\r", ((ix-1)*ny*nz+(iy-1)*nz+iz-1)/(nz*ny*nx)*100))
+            printfmt("{:.2f}\r", ((ix-1)*ny*nz+(iy-1)*nz+iz-1)/(nz*ny*nx)*100)
         end
         u_ana
     end
@@ -270,5 +269,4 @@ end
 
     close(M)
     @info "be patient, the garbage collection is probably deleting large scratch-files"
-end
-=#
\ No newline at end of file
+end
\ No newline at end of file
diff --git a/test/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl b/test/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl
index 4b6f83a..db5be8f 100644
--- a/test/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl
+++ b/test/jop_prop3DAcoVTIDenQ_DEO2_FDTD.jl
@@ -155,8 +155,7 @@ end
     close(F)
 end
 
-#=
-@testset "JopProp3DAcoVTIDenQ_DEO2_FDTD -- analytic, direct, for model type $modeltype, interpmethod=$interpmethod" for modeltype in ("vϵη","v"), interpmethod in (:linear,:hicks)
+@test_skip @testset "JopProp3DAcoVTIDenQ_DEO2_FDTD -- analytic, direct, for model type $modeltype, interpmethod=$interpmethod" for modeltype in ("vϵη","v"), interpmethod in (:linear,:hicks)
     z,y,x,dz,dy,dx,dtrec,dtmod,tmax,sz,sy,sx,c = 2000.0,500.0,2500.0,40.0,40.0,40.0,0.016,0.004,2.0,1000.0,250.0,20.0,1500.0
     nz,ny,nx,nt = round(Int,z/dz)+1,round(Int,y/dy)+1,round(Int,x/dx)+1,round(Int,tmax/dtrec)+1
 
@@ -187,7 +186,7 @@ end
                 end
             end
             u_ana[:,iz,iy,ix] = irfft(U,nt_pad)[1:nt]
-            printfmt("{:.2f}", ((ix-1)*ny*nz+(iy-1)*nz+iz-1)/(nz*ny*nx)*100))
+            printfmt("{:.2f}", ((ix-1)*ny*nz+(iy-1)*nz+iz-1)/(nz*ny*nx)*100)
         end
         u_ana
     end
@@ -268,5 +267,4 @@ end
 
     close(M)
     @info "be patient, the garbage collection is probably deleting large scratch-files"
-end
-=#
\ No newline at end of file
+end
\ No newline at end of file
diff --git a/test2D.jl b/test2D.jl
new file mode 100644
index 0000000..da58d0c
--- /dev/null
+++ b/test2D.jl
@@ -0,0 +1,175 @@
+using Revise
+using PyPlot
+using Formatting, FFTW, Jets, JetPackWaveFD, LinearAlgebra, SpecialFunctions, Statistics, Test, WaveFD
+
+modeltypes = (WaveFD.Prop2DAcoIsoDenQ_DEO2_FDTD_Model_V, WaveFD.Prop2DAcoIsoDenQ_DEO2_FDTD_Model_B, WaveFD.Prop2DAcoIsoDenQ_DEO2_FDTD_Model_VB)
+
+npad = 50;
+nx = 300;
+nz = 300;
+nt = 1501
+
+dx,dz,dt = 25.0,25.0,0.002
+xmin,zmin,tmin = 0.0,0.0,0.0
+xmax,zmax,tmax = xmin+dx*(nx-1),zmin+dz*(nz-1),tmin+dt*(nt-1)
+sx = dx * div(nx,2)
+sz = 1 * dz
+rx = dx*[npad:nx-npad-1;]
+rz = 2 * dz .* ones(length(rx))
+
+function make_op(interpmethod, modeltype, fs, v, b; comptype = Float32)
+    wavelet = WaveletCausalRicker(f=5.0)
+
+    local m
+    if modeltype == WaveFD.Prop2DAcoIsoDenQ_DEO2_FDTD_Model_V
+        kwargs = (b = reshape(b, nz, nx), )
+        m = reshape(v, nz, nx, 1)
+    elseif modeltype == WaveFD.Prop2DAcoIsoDenQ_DEO2_FDTD_Model_VB
+        kwargs = (nz = nz, nx = nx)
+        m = zeros(Float32, nz, nx, 2)
+        m[:,:,1] .= reshape(v, nz, nx)
+        m[:,:,2] .= reshape(b, nz, nx)
+    elseif modeltype == WaveFD.Prop2DAcoIsoDenQ_DEO2_FDTD_Model_B
+        kwargs = (v = reshape(v, nz, nx), )
+        m = reshape(b, nz, nx, 1)
+    end
+
+    F = JopNlProp2DAcoIsoDenQ_DEO2_FDTD(; kwargs..., 
+        z0 = zmin, x0 = xmin, dx = dx, dz = dz, sx = sx, sz = sz, rx = rx, rz = rz,
+        dtrec = dt, dtmod = dt, ntrec = nt, nsponge = npad, comptype = comptype, compscale = 1e-4,
+        freqQ = 5.0, qMin = 0.1, qInterior = 100.0, wavelet = wavelet, freesurface = fs, 
+        reportinterval = 0, interpmethod = interpmethod, isinterior = false)
+
+    m,F
+end
+
+boxsize = 10
+
+kz0 = div(nz,10)
+kz1 = kz0 - div(boxsize,2)
+kz2 = kz0 + div(boxsize,2)
+
+vx0 = 2 * div(nx,5)
+vx1 = vx0 - div(boxsize,2)
+vx2 = vx0 + div(boxsize,2)
+
+bx0 = 3 * div(nx,5)
+bx1 = bx0 - div(boxsize,2)
+bx2 = bx0 + div(boxsize,2)
+
+b1 =    1 .* ones(Float32, nz, nx, 1);
+b2 =    1 .* ones(Float32, nz, nx, 1); 
+v1 = 1500 .* ones(Float32, nz, nx, 1);
+v2 = 1500 .* ones(Float32, nz, nx, 1); 
+
+v2[kz1:kz2,vx1:vx2] .*= 1.25;
+b2[kz1:kz2,bx1:bx2] .*= 1.50;
+
+v3 = v2 .- v1
+b3 = b2 .- b1
+
+m1, F1 = make_op(:linear, WaveFD.Prop2DAcoIsoDenQ_DEO2_FDTD_Model_V, false, v1, b1);
+m2, F2 = make_op(:linear, WaveFD.Prop2DAcoIsoDenQ_DEO2_FDTD_Model_V, false, v2, b1);
+
+m3, F3 = make_op(:linear, WaveFD.Prop2DAcoIsoDenQ_DEO2_FDTD_Model_B, false, v1, b1);
+m4, F4 = make_op(:linear, WaveFD.Prop2DAcoIsoDenQ_DEO2_FDTD_Model_B, false, v1, b2);
+
+m5, F5 = make_op(:linear, WaveFD.Prop2DAcoIsoDenQ_DEO2_FDTD_Model_VB, false, v1, b1);
+m6, F6 = make_op(:linear, WaveFD.Prop2DAcoIsoDenQ_DEO2_FDTD_Model_VB, false, v2, b2);
+
+@show size(m1)
+@show size(m3)
+@show size(m5)
+
+@show extrema(m1[:,:,1]), extrema(m2[:,:,1])
+@show extrema(m3[:,:,1]), extrema(m4[:,:,1])
+@show extrema(m5[:,:,1]), extrema(m6[:,:,1])
+
+local d1,d2,d3,d4,d5,d6
+
+tv = @elapsed begin
+    @show "velocity only -- v1" 
+    d1 = F1 * m1;
+    @show "velocity only -- v2" 
+    d2 = F2 * m2;
+    @show "velocity only -- jacobian" 
+    J1 = jacobian!(F1,m1);
+    δd2 = J1 * (m2 .- m1);
+end
+@show tv
+
+tb = @elapsed begin
+    @show "buoyancy only -- b1" 
+    d3 = F3 * m3;
+    @show "buoyancy only -- b2" 
+    d4 = F4 * m4;
+    @show "buoyancy only -- jacobian" 
+    J3 = jacobian!(F3,m3);
+    δd3 = J3 * (m4 .- m3);
+end
+@show tb
+
+tvb = @elapsed begin
+    @show "velocity + buoyancy -- vb1" 
+    d5 = F5 * m5;
+    @show "velocity + buoyancy -- vb2" 
+    d6 = F6 * m6;
+    @show "velocity + buoyancy -- jacobian" 
+    J5 = jacobian!(F5,m5);
+    δd5 = J5 * (m6 .- m5);
+end
+@show tvb
+
+vmax = maximum(abs,v2); vmin = - vmax 
+bmax = maximum(abs,b2); bmin = - bmax 
+dmax = 0.25 * maximum(abs,d2); dmin = - dmax 
+
+close("all"); 
+
+shape = (nz, nx)
+
+figure(1,figsize=(18,12))
+subplot(2,3,1); imshow(reshape(v1, shape), cmap="seismic", aspect="auto", vmin=vmin, vmax=vmax); title("2D Velocity Background v1"); colorbar();
+subplot(2,3,2); imshow(reshape(v2, shape), cmap="seismic", aspect="auto", vmin=vmin, vmax=vmax); title("2D Velocity Perturbed v2"); colorbar();
+subplot(2,3,3); imshow(reshape(v3, shape), cmap="seismic", aspect="auto", vmin=vmin, vmax=vmax); title("2D Velocity Difference (v2-v1)"); colorbar();
+subplot(2,3,4); imshow(reshape(b1, shape), cmap="seismic", aspect="auto", vmin=bmin, vmax=bmax); title("2D Buoyancy Background b1"); colorbar();
+subplot(2,3,5); imshow(reshape(b2, shape), cmap="seismic", aspect="auto", vmin=bmin, vmax=bmax); title("2D Buoyancy Perturbed b2"); colorbar();
+subplot(2,3,6); imshow(reshape(b3, shape), cmap="seismic", aspect="auto", vmin=bmin, vmax=bmax); title("2D Buoyancy Difference (b2-b1)"); colorbar();
+tight_layout()
+display(gcf())
+savefig("figure.2D.model.png")
+
+scale = 10
+@show scale
+
+figure(3,figsize=(18,12))
+subplot(2,2,1); imshow(d1, cmap="seismic", aspect="auto", vmin=dmin, vmax=dmax); title("2D Reference F*v1"); colorbar();
+subplot(2,2,2); imshow(d2, cmap="seismic", aspect="auto", vmin=dmin, vmax=dmax); title("2D Velocity Only F*v2"); colorbar();
+subplot(2,2,3); imshow(d4, cmap="seismic", aspect="auto", vmin=dmin, vmax=dmax); title("2D Buoyancy Only F*b2"); colorbar();
+subplot(2,2,4); imshow(d6, cmap="seismic", aspect="auto", vmin=dmin, vmax=dmax); title("2D Velocity + Buoyancy F*vb2"); colorbar();
+tight_layout()
+display(gcf())
+savefig("figure.2D.data.nonlinear.png")
+
+figure(4,figsize=(18,12))
+subplot(2,2,2); imshow(scale .* (d2 .- d1), cmap="seismic", aspect="auto", vmin=dmin, vmax=dmax); title("2D Velocity Only F*v2 - F*v1 (($(scale)x)"); colorbar();
+subplot(2,2,3); imshow(scale .* (d4 .- d3), cmap="seismic", aspect="auto", vmin=dmin, vmax=dmax); title("2D Buoyancy Only F*b2 - F*b1 (($(scale)x)"); colorbar();
+subplot(2,2,4); imshow(scale .* (d6 .- d5), cmap="seismic", aspect="auto", vmin=dmin, vmax=dmax); title("2D Velocity + Buoyancy F*vb2 - F*vb1 (($(scale)x)"); colorbar();
+tight_layout()
+display(gcf())
+savefig("figure.2D.data.difference.png")
+
+figure(5,figsize=(18,12))
+subplot(2,2,2); imshow(scale .* (δd2), cmap="seismic", aspect="auto", vmin=dmin, vmax=dmax); title("2D Velocity Only J*(v2-v1) ($(scale)x)"); colorbar();
+subplot(2,2,3); imshow(scale .* (δd3), cmap="seismic", aspect="auto", vmin=dmin, vmax=dmax); title("2D Buoyancy Only J*(b2-b1) (($(scale)x)"); colorbar();
+subplot(2,2,4); imshow(scale .* (δd5), cmap="seismic", aspect="auto", vmin=dmin, vmax=dmax); title("2D Velocity + Buoyancy J*(vb2-vb1) (($(scale)x)"); colorbar();
+tight_layout()
+display(gcf())
+savefig("figure.2D.data.linear.png")
+
+close(F1)
+close(F2)
+close(F3)
+close(F4)
+close(F5)
+close(F6)
diff --git a/test3D.jl b/test3D.jl
new file mode 100644
index 0000000..14b8a34
--- /dev/null
+++ b/test3D.jl
@@ -0,0 +1,183 @@
+using Revise
+using PyPlot
+using Formatting, FFTW, Jets, JetPackWaveFD, LinearAlgebra, SpecialFunctions, Statistics, Test, WaveFD
+
+modeltypes = (WaveFD.Prop2DAcoIsoDenQ_DEO2_FDTD_Model_V, WaveFD.Prop2DAcoIsoDenQ_DEO2_FDTD_Model_B, WaveFD.Prop2DAcoIsoDenQ_DEO2_FDTD_Model_VB)
+
+npad = 20;
+nx = 220;
+ny = 210;
+nz = 200;
+nt = 1001
+
+dx,dy,dz,dt = 25.0,25.0,25.0,0.002
+xmin,ymin,zmin,tmin = 0.0,0.0,0.0,0.0
+xmax,ymax,zmax,tmax = xmin+dx*(nx-1),ymin+dy*(ny-1),zmin+dz*(nz-1),tmin+dt*(nt-1)
+sx = dx * div(nx,2)
+sy = dy * div(ny,2)
+sz = 1 * dz
+rx = dx*[npad:nx-npad-1;]
+ry = div(ny,2) * dy .* ones(length(rx))
+rz =        2  * dz .* ones(length(rx))
+
+function make_op(interpmethod, modeltype, fs, v, b; comptype = Float32)
+    wavelet = WaveletCausalRicker(f=5.0)
+
+    local m
+    if modeltype == WaveFD.Prop2DAcoIsoDenQ_DEO2_FDTD_Model_V
+        kwargs = (b = reshape(b, nz, ny, nx), )
+        m = reshape(v, nz, ny, nx, 1)
+    elseif modeltype == WaveFD.Prop2DAcoIsoDenQ_DEO2_FDTD_Model_VB
+        kwargs = (nz = nz, ny = ny, nx = nx)
+        m = zeros(Float32, nz, ny, nx, 2)
+        m[:,:,:,1] .= reshape(v,nz,ny,nx)
+        m[:,:,:,2] .= reshape(b,nz,ny,nx)
+    elseif modeltype == WaveFD.Prop2DAcoIsoDenQ_DEO2_FDTD_Model_B
+        kwargs = (v = reshape(v, nz, ny, nx), )
+        m = reshape(b, nz, ny, nx, 1)
+    end
+
+    F = JopNlProp3DAcoIsoDenQ_DEO2_FDTD(; kwargs..., 
+        z0 = zmin, y0 = ymin, x0 = xmin, dz = dz, dy = dy, dx = dx, 
+        sz = sz, sy = sy, sx = sx, rz = rz, ry = ry, rx = rx, 
+        dtrec = dt, dtmod = dt, ntrec = nt, nsponge = npad, comptype = comptype, compscale = 1e-4,
+        freqQ = 5.0, qMin = 0.1, qInterior = 100.0, wavelet = wavelet, freesurface = fs, 
+        reportinterval = 0, interpmethod = interpmethod, isinterior = false)
+
+    m,F
+end
+
+boxsize = 10
+
+kz0 = div(nz,10)
+kz1 = kz0 - div(boxsize,2)
+kz2 = kz0 + div(boxsize,2)
+
+ky0 = div(ny,2)
+ky1 = ky0 - div(boxsize,2)
+ky2 = ky0 + div(boxsize,2)
+
+vx0 = 2 * div(nx,5)
+vx1 = vx0 - div(boxsize,2)
+vx2 = vx0 + div(boxsize,2)
+
+bx0 = 3 * div(nx,5)
+bx1 = bx0 - div(boxsize,2)
+bx2 = bx0 + div(boxsize,2)
+
+b1 =    1 .* ones(Float32, nz, ny, nx, 1);
+b2 =    1 .* ones(Float32, nz, ny, nx, 1); 
+v1 = 1500 .* ones(Float32, nz, ny, nx, 1);
+v2 = 1500 .* ones(Float32, nz, ny, nx, 1); 
+
+v2[kz1:kz2,ky1:ky2,vx1:vx2] .*= 1.25;
+b2[kz1:kz2,ky1:ky2,bx1:bx2] .*= 1.50;
+
+v3 = v2 .- v1
+b3 = b2 .- b1
+
+m1, F1 = make_op(:linear, WaveFD.Prop2DAcoIsoDenQ_DEO2_FDTD_Model_V, false, v1, b1);
+m2, F2 = make_op(:linear, WaveFD.Prop2DAcoIsoDenQ_DEO2_FDTD_Model_V, false, v2, b1);
+
+m3, F3 = make_op(:linear, WaveFD.Prop2DAcoIsoDenQ_DEO2_FDTD_Model_B, false, v1, b1);
+m4, F4 = make_op(:linear, WaveFD.Prop2DAcoIsoDenQ_DEO2_FDTD_Model_B, false, v1, b2);
+
+m5, F5 = make_op(:linear, WaveFD.Prop2DAcoIsoDenQ_DEO2_FDTD_Model_VB, false, v1, b1);
+m6, F6 = make_op(:linear, WaveFD.Prop2DAcoIsoDenQ_DEO2_FDTD_Model_VB, false, v2, b2);
+
+@show size(m1)
+@show size(m3)
+@show size(m5)
+
+@show extrema(m1[:,:,:,1]), extrema(m2[:,:,:,1])
+@show extrema(m3[:,:,:,1]), extrema(m4[:,:,:,1])
+@show extrema(m5[:,:,:,1]), extrema(m6[:,:,:,1])
+
+local d1,d2,d3,d4,d5,d6
+
+tv = @elapsed begin
+    @show "velocity only -- v1" 
+    d1 = F1 * m1;
+    @show "velocity only -- v2" 
+    d2 = F2 * m2;
+    @show "velocity only -- jacobian" 
+    J1 = jacobian!(F1,m1);
+    δd2 = J1 * (m2 .- m1);
+end
+@show tv
+
+tb = @elapsed begin
+    @show "buoyancy only -- b1" 
+    d3 = F3 * m3;
+    @show "buoyancy only -- b2" 
+    d4 = F4 * m4;
+    @show "buoyancy only -- jacobian" 
+    J3 = jacobian!(F3,m3);
+    δd3 = J3 * (m4 .- m3);
+end
+@show tb
+
+tvb = @elapsed begin
+    @show "velocity + buoyancy -- vb1" 
+    d5 = F5 * m5;
+    @show "velocity + buoyancy -- vb2" 
+    d6 = F6 * m6;
+    @show "velocity + buoyancy -- jacobian" 
+    J5 = jacobian!(F5,m5);
+    δd5 = J5 * (m6 .- m5);
+end
+@show tvb
+
+vmax = maximum(abs,v2); vmin = - vmax 
+bmax = maximum(abs,b2); bmin = - bmax 
+dmax = 0.05 * maximum(abs,d2); dmin = - dmax 
+
+close("all"); 
+
+shape = (nz, nx)
+
+figure(1,figsize=(18,12))
+subplot(2,3,1); imshow(reshape(v1[:,div(ny,2),:], shape), cmap="seismic", aspect="auto", vmin=vmin, vmax=vmax); title("3D Velocity Background v1"); colorbar();
+subplot(2,3,2); imshow(reshape(v2[:,div(ny,2),:], shape), cmap="seismic", aspect="auto", vmin=vmin, vmax=vmax); title("3D Velocity Perturbed v2"); colorbar();
+subplot(2,3,3); imshow(reshape(v3[:,div(ny,2),:], shape), cmap="seismic", aspect="auto", vmin=vmin, vmax=vmax); title("3D Velocity Difference (v2-v1)"); colorbar();
+subplot(2,3,4); imshow(reshape(b1[:,div(ny,2),:], shape), cmap="seismic", aspect="auto", vmin=bmin, vmax=bmax); title("3D Buoyancy Background b1"); colorbar();
+subplot(2,3,5); imshow(reshape(b2[:,div(ny,2),:], shape), cmap="seismic", aspect="auto", vmin=bmin, vmax=bmax); title("3D Buoyancy Perturbed b2"); colorbar();
+subplot(2,3,6); imshow(reshape(b3[:,div(ny,2),:], shape), cmap="seismic", aspect="auto", vmin=bmin, vmax=bmax); title("3D Buoyancy Difference (b2-v1)"); colorbar();
+tight_layout()
+display(gcf())
+savefig("figure.3D.model.png")
+
+scale = 10
+@show scale
+
+figure(3,figsize=(18,12))
+subplot(2,2,1); imshow(d1, cmap="seismic", aspect="auto", vmin=dmin, vmax=dmax); title("3D Reference F*v1"); colorbar();
+subplot(2,2,2); imshow(d2, cmap="seismic", aspect="auto", vmin=dmin, vmax=dmax); title("3D Velocity Only F*v2"); colorbar();
+subplot(2,2,3); imshow(d4, cmap="seismic", aspect="auto", vmin=dmin, vmax=dmax); title("3D Buoyancy Only F*b2"); colorbar();
+subplot(2,2,4); imshow(d6, cmap="seismic", aspect="auto", vmin=dmin, vmax=dmax); title("3D Velocity + Buoyancy F*vb2"); colorbar();
+tight_layout()
+display(gcf())
+savefig("figure.3D.data.nonlinear.png")
+
+figure(4,figsize=(18,12))
+subplot(2,2,2); imshow(scale .* (d2 .- d1), cmap="seismic", aspect="auto", vmin=dmin, vmax=dmax); title("3D Velocity Only F*v2 - F*v1 (($(scale)x)"); colorbar();
+subplot(2,2,3); imshow(scale .* (d4 .- d3), cmap="seismic", aspect="auto", vmin=dmin, vmax=dmax); title("3D Buoyancy Only F*b2 - F*b1 (($(scale)x)"); colorbar();
+subplot(2,2,4); imshow(scale .* (d6 .- d5), cmap="seismic", aspect="auto", vmin=dmin, vmax=dmax); title("3D Velocity + Buoyancy F*vb2 - F*vb1 (($(scale)x)"); colorbar();
+tight_layout()
+display(gcf())
+savefig("figure.3D.data.difference.png")
+
+figure(5,figsize=(18,12))
+subplot(2,2,2); imshow(scale .* (δd2), cmap="seismic", aspect="auto", vmin=dmin, vmax=dmax); title("3D Velocity Only J*(v2-v1) ($(scale)x)"); colorbar();
+subplot(2,2,3); imshow(scale .* (δd3), cmap="seismic", aspect="auto", vmin=dmin, vmax=dmax); title("3D Buoyancy Only J*(b2-b1) (($(scale)x)"); colorbar();
+subplot(2,2,4); imshow(scale .* (δd5), cmap="seismic", aspect="auto", vmin=dmin, vmax=dmax); title("3D Velocity + Buoyancy J*(vb2-vb1) (($(scale)x)"); colorbar();
+tight_layout()
+display(gcf())
+savefig("figure.3D.data.linear.png")
+
+close(F1)
+close(F2)
+close(F3)
+close(F4)
+close(F5)
+close(F6)