diff --git a/.github/workflows/ci-examples.yml b/.github/workflows/ci-examples.yml
index 10f5f1889..42e72a950 100644
--- a/.github/workflows/ci-examples.yml
+++ b/.github/workflows/ci-examples.yml
@@ -35,7 +35,7 @@ jobs:
   run-examples:
     runs-on: ubuntu-latest
     needs: list-examples
-    name: ${{ matrix.example }} on Julia ${{ matrix.version }}
+    name: ${{ matrix.example }}
 
     env:
       DEVITO_ARCH: gcc-9
@@ -48,15 +48,8 @@ jobs:
       fail-fast: false
       matrix:
         example: ${{ fromJson(needs.list-examples.outputs.matrix) }}
-        version: ['1.8']
-
-        include:
-         - example: "examples/scripts/modeling_basic_2D.jl"
-           version: '1.6'
-
-         - example: "examples/scripts/modeling_basic_2D.jl"
-           version: '1.7'
-
+        version: ['1']
+  
     steps:
       - name: Checkout JUDI
         uses: actions/checkout@v3
diff --git a/.github/workflows/ci-judi.yml b/.github/workflows/ci-judi.yml
index 3fc1c9221..1e97c04f9 100644
--- a/.github/workflows/ci-judi.yml
+++ b/.github/workflows/ci-judi.yml
@@ -28,7 +28,7 @@ jobs:
       fail-fast: false
 
       matrix:
-        version: ['1.6', '1.7', '1.8', '1.9']
+        version: ['1.6', '1.7', '1.8', '1.9', '1.10']
         os: [ubuntu-latest, macos-latest]
 
     steps:
diff --git a/.github/workflows/ci-op.yml b/.github/workflows/ci-op.yml
index e1ef93fee..a18d7a443 100644
--- a/.github/workflows/ci-op.yml
+++ b/.github/workflows/ci-op.yml
@@ -37,27 +37,27 @@ jobs:
         include:
           - os: macos-latest
             version: '1.6'
-            op: "BASICS"
+            op: "ISO_OP"
             omp: 1
 
           - os: macos-latest
-            version: '1.7'
-            op: "BASICS"
+            version: '1.8'
+            op: "ISO_OP"
             omp: 1
 
           - os: macos-latest
-            version: '1.8'
-            op: "BASICS"
+            version: '1.9'
+            op: "ISO_OP"
             omp: 1
-          
+
           - os: ubuntu-latest
             version: '1.9'
             op: "VISCO_AC_OP"
             omp: 2
 
           - os: ubuntu-latest
-            version: '1.6'
-            op: "BASICS"
+            version: '1.10'
+            op: "ISO_OP"
             omp: 2
 
     steps:
diff --git a/Project.toml b/Project.toml
index e0490035c..ce6ef2777 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "JUDI"
 uuid = "f3b833dc-6b2e-5b9c-b940-873ed6319979"
 authors = ["Philipp Witte, Mathias Louboutin"]
-version = "3.3.8"
+version = "3.3.9"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
diff --git a/examples/scripts/lsrtm_2D.jl b/examples/scripts/lsrtm_2D.jl
index 0d6ab6272..d4f41ce67 100644
--- a/examples/scripts/lsrtm_2D.jl
+++ b/examples/scripts/lsrtm_2D.jl
@@ -12,15 +12,19 @@ if ~isfile("$(JUDI.JUDI_DATA)/marmousi_model.h5")
 end
 n, d, o, m0 = read(h5open("$(JUDI.JUDI_DATA)/marmousi_model.h5", "r"), "n", "d", "o", "m0")
 
+# Smaller model for CI
+if get(ENV, "NITER", "10") == "2"
+    fact = 4
+    m0 = m0[1:fact:end, 1:fact:end]
+    n = size(m0)
+    d = d .* fact
+else
+    fact = 1
+end
+
 # Set up model structure
 model0 = Model(n, d, o, m0)
-grad_mem = 40 # Based on n and CFL condition
-
-# Coarsen for CI
-if get(ENV, "CI", nothing) == "true"
-    model0 = Model(ceil.(Int, n ./ 2), d .* 2, o, m0[1:2:end, 1:2:end])
-    grad_mem = 5
-end
+grad_mem = 40 / (fact^3) # Based on n and CFL condition
 
 # Load data
 if ~isfile("$(JUDI.JUDI_DATA)/marmousi_2D.segy")
@@ -37,7 +41,7 @@ q = judiVector(src_geometry, wavelet)
 ###################################################################################################
 # Infer subsampling based on free memory
 mem = Sys.free_memory()/(1024^3)
-t_sub = max(1, ceil(Int, nworkers()*grad_mem/mem))
+t_sub = max(1, ceil(Int, .5*nworkers()*grad_mem/mem))
 
 # Setup operators
 opt = Options(subsampling_factor=t_sub, isic=true)  # ~40 GB of memory per source without subsampling
@@ -53,7 +57,7 @@ Ml = judiDataMute(q.geometry, d_lin.geometry)
 #' set up number of iterations
 niter = parse(Int, get(ENV, "NITER", "10"))
 # Default to 64, 5 for CI only with NITER=1
-nsrc = 5 * parse(Int, get(ENV, "NITER", "$(q.nsrc ÷ 5)"))
+nsrc = 5 * parse(Int, get(ENV, "NITER", "10"))
 indsrc = randperm(q.nsrc)[1:nsrc]
 lsqr_sol = zeros(Float32, prod(model0.n))
 
diff --git a/examples/scripts/splsrtm_2D.jl b/examples/scripts/splsrtm_2D.jl
index dea1c21bd..46c80cb61 100644
--- a/examples/scripts/splsrtm_2D.jl
+++ b/examples/scripts/splsrtm_2D.jl
@@ -12,9 +12,19 @@ if ~isfile("$(JUDI.JUDI_DATA)/marmousi_model.h5")
 end
 n, d, o, m0 = read(h5open("$(JUDI.JUDI_DATA)/marmousi_model.h5", "r"), "n", "d", "o", "m0")
 
+# Smaller model for CI
+if get(ENV, "NITER", "10") == "2"
+    fact = 4
+    m0 = m0[1:fact:end, 1:fact:end]
+    n = size(m0)
+    d = d .* fact
+else
+    fact = 1
+end
+
 # Set up model structure
 model0 = Model(n, d, o, m0)
-grad_mem = 40 # Based on n and CFL condition
+grad_mem = 40 / (fact^3) # Based on n and CFL condition
 
 # Coarsen for CI
 if get(ENV, "CI", nothing) == "true"
@@ -37,7 +47,7 @@ q = judiVector(src_geometry, wavelet)
 ###################################################################################################
 # Infer subsampling based on free memory
 mem = Sys.free_memory()/(1024^3)
-t_sub = max(1, ceil(Int, nworkers()*grad_mem/mem))
+t_sub = max(1, ceil(Int, .5*nworkers()*grad_mem/mem))
 
 # Setup operators
 opt = Options(subsampling_factor=t_sub, isic=true)  # ~40 GB of memory per source without subsampling
@@ -54,7 +64,7 @@ C = joEye(prod(model0.n); DDT=Float32, RDT=Float32)
 # If available use curvelet instead for better result
 
 # Setup linearized bregman
-batchsize = 5 * parse(Int, get(ENV, "NITER", "$(q.nsrc ÷ 5)"))
+batchsize = 5 * parse(Int, get(ENV, "NITER", "10"))
 niter = parse(Int, get(ENV, "NITER", "10"))
 g_scale = 0
 
@@ -65,7 +75,7 @@ function obj(x)
 
     residual = Ml[inds]*J[inds]*Mr*dm - Ml[inds]*d_lin[inds]
     # grad
-    G = reshape(Mr'J[inds]'*Ml[inds]'*residual, model0.n)
+    G = reshape(Mr'*J[inds]'*Ml[inds]'*residual, model0.n)
     g_scale == 0 && (global g_scale = .05f0/maximum(G))
     G .*= g_scale
     return .5f0*norm(residual)^2, G[:]
diff --git a/src/JUDI.jl b/src/JUDI.jl
index 8889aa489..4365b07d1 100644
--- a/src/JUDI.jl
+++ b/src/JUDI.jl
@@ -24,7 +24,7 @@ import Base.depwarn
 import Base.*, Base./, Base.+, Base.-, Base.==, Base.\
 import Base.copy!, Base.copy, Base.copyto!, Base.deepcopy, Base.summary
 import Base.sum, Base.ndims, Base.reshape, Base.fill!, Base.axes, Base.dotview
-import Base.eltype, Base.length, Base.size, Base.iterate, Base.show, Base.display, Base.showarg
+import Base.eltype, Base.length, Base.size, Base.iterate, Base.show, Base.display, Base.showarg, Base.elsize
 import Base.maximum, Base.minimum, Base.push!
 import Base.Broadcast.ArrayStyle, Base.Broadcast.extrude
 import Base.Broadcast.broadcasted, Base.BroadcastStyle, Base.Broadcast.DefaultArrayStyle, Base.Broadcast, Base.broadcast!
diff --git a/src/TimeModeling/Preconditioners/ModelPreconditioners.jl b/src/TimeModeling/Preconditioners/ModelPreconditioners.jl
index 0e617306a..bf0e9cf5d 100644
--- a/src/TimeModeling/Preconditioners/ModelPreconditioners.jl
+++ b/src/TimeModeling/Preconditioners/ModelPreconditioners.jl
@@ -69,7 +69,10 @@ function matvec(D::TopMute{T, N}, x::Array{T, N}) where {T, N}
     taper = D.taperwidth < 2 ? 1 : _taper(Val(:reflection), D.taperwidth)
     for i in CartesianIndices(D.wb)
         out[i, 1:D.wb[i]-D.taperwidth] .= 0
-        out[i, D.wb[i]-D.taperwidth+1:D.wb[i]] .*= taper
+        s = max(D.wb[i]-D.taperwidth+1, 1)
+        ni = length(s:D.wb[i])
+        tap = D.taperwidth < 2 ? 1 : taper[end-ni+1:end]
+        out[i, s:D.wb[i]] .*= tap
     end
     out
 end
@@ -182,7 +185,7 @@ end
 function set_val(I::judiIllumination{T, M, K, R}, mode, v) where {T, M, K, R}
     key = mode ∈ [:forward, :born] ? "u" : "v"
     if key in keys(I.illums)
-        I.illums[key] .= v
+        combine!(I.illums[key], v)
     end
 end
 
diff --git a/src/TimeModeling/Types/ModelStructure.jl b/src/TimeModeling/Types/ModelStructure.jl
index 68e3e1af3..66819c29b 100644
--- a/src/TimeModeling/Types/ModelStructure.jl
+++ b/src/TimeModeling/Types/ModelStructure.jl
@@ -144,6 +144,8 @@ dotview(m::PhysicalParameter, i) = Base.dotview(m.data, i)
 getindex(A::PhysicalParameter{T, N}, i::Int) where {T<:Real, N} = A.data[i]
 getindex(A::PhysicalParameter{T, N}, ::Colon) where {T<:Real, N} = A.data[:]
 
+elsize(A::PhysicalParameter)  = elsize(A.data)
+
 get_step(r::StepRange) = r.step
 get_step(r) = 1
 
@@ -192,12 +194,20 @@ function combine(op, A::PhysicalParameter{T, N}, B::PhysicalParameter{T, N}) whe
     mn = max.(ea, eb)
     ia = [s:e for (s, e) in zip(sa, ea)]
     ib = [s:e for (s, e) in zip(sb, eb)]
-    out = zeros(T, mn)
-    out[ia...] .= A.data
-    broadcast!(op, view(out, ib...),  view(out, ib...), B.data)
-    return PhysicalParameter(mn, A.d, o, out)
+    if isnothing(op)
+        @assert A.n == mn
+        A.data[ib...] .= B.data
+        return nothing
+    else
+        out = zeros(T, mn)
+        out[ia...] .= A.data
+        broadcast!(op, view(out, ib...),  view(out, ib...), B.data)
+        return PhysicalParameter(mn, A.d, o, out)
+    end
 end
 
+combine!(A::PhysicalParameter{T, N}, B::PhysicalParameter{T, N}) where {T<:Real, N} = combine(nothing, A, B)
+
 for op in [:+, :-, :*, :/]
     @eval function $(op)(A::PhysicalParameter{T, N}, B::PhysicalParameter{T, N}) where {T<:Real, N}
         if compare(A, B)
diff --git a/src/TimeModeling/Types/broadcasting.jl b/src/TimeModeling/Types/broadcasting.jl
index 5a18d93e3..72df0a55e 100644
--- a/src/TimeModeling/Types/broadcasting.jl
+++ b/src/TimeModeling/Types/broadcasting.jl
@@ -27,7 +27,7 @@ end
 
 # Broadcasted custom type
 check_compat() = true
-get_src(ms::judiMultiSourceVector, j) = ms.data[j]
+get_src(ms::judiMultiSourceVector, j) = make_input(ms[j])
 get_src(v::Vector{<:Array}, j) = v[j]
 get_src(v::Array{T, N}, j) where {T<:Number, N} = v
 get_src(n::Number, j) = n
diff --git a/src/pysource/fields.py b/src/pysource/fields.py
index be3a6e4ef..55953f7ac 100644
--- a/src/pysource/fields.py
+++ b/src/pysource/fields.py
@@ -136,7 +136,7 @@ def wavefield_subsampled(model, u, nt, t_sub, space_order=8):
     for wf in as_tuple(u):
         usave = TimeFunction(name='us_%s' % wf.name, grid=model.grid, time_order=2,
                              space_order=space_order, time_dim=time_subsampled,
-                             save=nsave)
+                             save=int(nsave))
         wf_s.append(usave)
     return wf_s
 
diff --git a/src/pysource/propagators.py b/src/pysource/propagators.py
index 9d6042522..109b6b760 100644
--- a/src/pysource/propagators.py
+++ b/src/pysource/propagators.py
@@ -97,6 +97,8 @@ def gradient(model, residual, rcv_coords, u, return_op=False, space_order=8, fw=
     v = wavefield(model, space_order, fw=not fw)
     try:
         t_sub = as_tuple(u)[0].indices[0]._factor
+        if isinstance(t_sub, Constant):
+            t_sub = t_sub.data
     except AttributeError:
         t_sub = 1