diff --git a/scripts/mit_inv.jl b/scripts/mit_inv.jl
index 633a6315..4d3b805c 100644
--- a/scripts/mit_inv.jl
+++ b/scripts/mit_inv.jl
@@ -40,12 +40,13 @@ flush(stdout)
 # The first guess for the tracer concentration should be close to the actual tracer concentration
 # take first guess as θtrue
 cvec=(vec(θtrue))
-θguess = unvec(θtrue,cvec)
-
+#θguess = unvec(θtrue,cvec)
 
+ucnew = load("ucnew.jld")["ucnew"]
 #first guess tracer control vector is near zero, and we want this to remain relatively small
-u =  Field(-0.01.*ones(size(γ.wet)),γ,θtrue.name,θtrue.longname,θtrue.units)
-uvec = vec(u)
+u1 =  Field(0.0.*ones(size(γ.wet)),γ,θtrue.name,θtrue.longname,θtrue.units)
+uvec = ucnew
+u =  unvec(u1,ucnew)
 q = zeros(size(uvec))
 q[surfind] = cvec[surfind]
 
@@ -63,8 +64,10 @@ Q⁻ = Diagonal(Qdiag)
 
 ind_1 = load("ind_1.jld")["ind_1"]
 ind_2 = load("ind_2.jld")["ind_2"]
+ufnew = load("ufnew.jld")["ufnew"]
+
 
-A0 = sparse(ind_1, ind_2, ones(length(ind_1))*0.1)
+Adummy = sparse(ind_1, ind_2, ones(length(ind_1))*0.1)
 
 
 #test: apply A0 a bunch of times
@@ -91,9 +94,11 @@ A0 = sparse(ind_1, ind_2, ones(length(ind_1))*0.1)
 #readline()
 
 
-non_zero_indices1, non_zero_indices2, non_zero_values = findnz(A0)
+non_zero_indices1, non_zero_indices2, non_zero_values = findnz(Adummy)
 
+Afguess = sparse(non_zero_indices1, non_zero_indices2, ufnew)
 
+non_zero_indices1, non_zero_indices2, non_zero_values = findnz(Afguess)
 
 non_zero_indices = hcat(non_zero_indices1, non_zero_indices2)
 
@@ -105,9 +110,9 @@ print(ulength)
 print(sum(γ.wet))
 
 # get sample J value
-F = costfunction_gridded_model(convec,non_zero_indices,u,A0,ctrue,cvec,q,W⁻,Q⁻,γ)
-fg!(F,G,x) = costfunction_gridded_model!(F,G,x,non_zero_indices,u,A0,ctrue,cvec,q,W⁻,Q⁻,γ)
-fg(x) = costfunction_gridded_model(x,non_zero_indices,u,A0,ctrue,cvec,q,W⁻,Q⁻,γ)
+F = costfunction_gridded_model(convec,non_zero_indices,u,Adummy,ctrue,cvec,q,W⁻,Q⁻,γ)
+fg!(F,G,x) = costfunction_gridded_model!(F,G,x,non_zero_indices,u,Adummy,ctrue,cvec,q,W⁻,Q⁻,γ)
+fg(x) = costfunction_gridded_model(x,non_zero_indices,u,Adummy,ctrue,cvec,q,W⁻,Q⁻,γ)
 f(x) = fg(x)[1]
 J₀,gJ₀ = fg(convec)
 
@@ -115,8 +120,9 @@ J₀,gJ₀ = fg(convec)
 # check with forward differences
 ϵ = 1e-3
 println(size(length(convec)))
-ii = rand(1:length(convec))
-ii=748529+10000
+#ii = rand(1:length(convec))
+#ii=748529+10000
+ii=4000
 println("Location for test =",ii)
 flush(stdout)
 δu = copy(convec); δu[ii] += ϵ
@@ -132,24 +138,24 @@ println("Percent error ",100*abs(∇f - ∇f_finite)/abs(∇f + ∇f_finite))
 flush(stdout)
 #### end gradient check #################
 
-#print(length(convec))
+print(length(convec))
 # filter the data with an Optim.jl method
-iterations = 2
+iterations = 20000
 out = steadyclimatology_optim(convec,fg!,iterations)
 
 # reconstruct by hand to double-check.
-ũ = unvec(θguess,Float32.(out.minimizer[begin:ulength]))
+#ũ = unvec(θguess,Float32.(out.minimizer[begin:ulength]))
 
 #println(typeof(ũ))
 
 # reconstruct tracer map
-c₀ = θguess
-c̃  = θguess+ũ
+#c₀ = θguess
+#c̃  = θguess+ũ
 
-Δc̃ = c̃ - θtrue
-Δc₀ = θguess - θtrue
+#Δc̃ = c̃ - θtrue
+#Δc₀ = θguess - θtrue
 
-Anew = A0 + sparse(non_zero_indices[:, 1], non_zero_indices[:, 2], out.minimizer[ulength+1:end])
+Anew = Adummy + sparse(non_zero_indices[:, 1], non_zero_indices[:, 2], out.minimizer[ulength+1:end])
 
 maxA = maximum(Anew)
 minA = minimum(Anew)
@@ -157,14 +163,23 @@ minA = minimum(Anew)
 println("Max A:$maxA")
 println("Min A:$minA")
 
+ucout = out.minimizer[begin:ulength]
+
+maxuc = maximum(ucout)
+minuc = minimum(ucout)
+
+println("Max uc:$maxuc")
+println("Min uc:$minuc")
+
+
 flush(stdout)
 onesvec = ones(size(q))
 
 #oldf = sum((non_zero_values).^2)
 #newf = sum((out.minimizer[ulength+1:end]).^2)
-tracer_cons1 = sum((Q⁻*A0*cvec - q).^2)
-tracer_cons2 = sum((Q⁻*Anew*(cvec+out.minimizer[begin:ulength]) - q).^2)
-mass_cons1 = sum((A0*onesvec).^2)
+tracer_cons1 = sum(((Adummy+Afguess)*(cvec+u) - q).^2)
+tracer_cons2 = sum((Anew*(cvec+out.minimizer[begin:ulength]) - q).^2)
+mass_cons1 = sum(((Adummy+ Afguess)*onesvec).^2)
 mass_cons2 = sum((Anew*onesvec).^2)
 
 println("old tracer cons:$tracer_cons1")
@@ -172,7 +187,14 @@ println("new tracer cons:$tracer_cons2")
 println("old mass cons:$mass_cons1")
 println("new mass cons:$mass_cons2")
 
-#b = getsurfaceboundary(θtrue)
+b = getsurfaceboundary(θtrue)
+
+save("ufnew.jld","ufnew", out.minimizer[ulength+1:end])
+save("ucnew.jld","ucnew", out.minimizer[begin:ulength])
+
+
+
+
 
 #A0
 
@@ -181,7 +203,22 @@ println("new mass cons:$mass_cons2")
 #Alu2 = lu(Anew)
 
 
-#final_θ = steadyinversion(Anew,b,γ)
+final_θ = steadyinversion(Anew,b,γ)
+#
+
+c1 = (vec(final_θ))
+
+test = sum(Anew * c1 - q)
+
+
+max_θ = maximum(filter(!isnan,final_θ.tracer))
+min_θ = minimum(filter(!isnan,final_θ.tracer))
+
+println("final_θ max:$max_θ")
+println("final_θ min:$min_θ")
+println("test:$test")
+
+#writefield("final_theta.nc",final_θ)
 
 # plot the difference
 #level = 15 # your choice 1-33
diff --git a/src/TMI.jl b/src/TMI.jl
index d79ab05e..7d5e307b 100644
--- a/src/TMI.jl
+++ b/src/TMI.jl
@@ -1672,22 +1672,26 @@ function costfunction_gridded_model(convec,non_zero_indices,u₀::Field{T},A0,y:
     #dummy,dummy,fguess  = findnz(A0)
     #dummy,dummy,fnow  = findnz(A) 
     onesvec = ones(size(uvec))
-    #csum = Wⁱ * uvec+c
+    csum = uvec+c
 
 
     # find lagrange multipliers
-    #muk = transpose(A) * Qⁱ * (A * c - q)
-    dAcdf = spzeros(length(ufvec),length(c))
-    dA1df = spzeros(length(ufvec),length(c))
+    muk = transpose(A) * Qⁱ * (A * csum - q)
+    dAcdf = spzeros(length(ufvec))#,length(c))
+    dA1df = spzeros(length(ufvec))#,length(c))
+    A1 = (A * onesvec)
+    Acq = Qⁱ*(A * csum - q)
     for ii in eachindex(ufvec)
-          dAcdf[ii,non_zero_indices[ii, 1]] = c[non_zero_indices[ii, 2]]
-          dA1df[ii,non_zero_indices[ii, 1]] = 1
+          #dAcdf[ii,non_zero_indices[ii, 1]] = c[non_zero_indices[ii, 2]]
+          #dA1df[ii,non_zero_indices[ii, 1]] = 1
+	  dAcdf[ii] = sum(csum[non_zero_indices[ii, 2]].*Acq[non_zero_indices[ii, 1]])
+	  dA1df[ii] = sum(A1[non_zero_indices[ii, 1]])
     end
     
 
-    dAdf_terms = dA1df * (A * onesvec) + dAcdf * Qⁱ * (A * c - q)# + dA1df * (A * onesvec)
+    dAdf_terms = dA1df + dAcdf #* Qⁱ * (A * c - q)# + dA1df * (A * onesvec)
 
-    J =  10^3 * uvec ⋅ uvec + transpose(A * c - q) * Qⁱ * (A*c - q) + 
+    J =  10 *uvec ⋅ uvec + transpose(A * csum - q) * Qⁱ * (A*csum - q) + 
     #2 * transpose(muk)*( Wⁱ * uvec+c-y)# + sum(ufvec)#sum(A*onesvec)
           transpose(A * onesvec)* (A* onesvec)# + 10 * ufvec ⋅ ufvec.^3
 
@@ -1698,7 +1702,7 @@ function costfunction_gridded_model(convec,non_zero_indices,u₀::Field{T},A0,y:
         if ii <= ulength 
           #this is the derivative of the cost function wrt the part of the control vector
           # associated with the tracer concentration
-          guvec[ii] =  2 *  10^3 * uvec[ii] #- (2 * transpose(muk) * Wⁱ)[ii]
+          guvec[ii] =  20 *  uvec[ii] + 2 * transpose(muk)[ii]
         else
           #this is the derivative of the cost function wrt the part of the control vector
           # associated with the transport vector
@@ -1723,23 +1727,27 @@ function costfunction_gridded_model!(J,guvec,convec::Vector{T},non_zero_indices,
     #dummy,dummy,fguess  = findnz(A0)
     #dummy,dummy,fnow  = findnz(A)
     onesvec = ones(size(uvec))
-    #csum = Wⁱ * uvec+c
+    csum = uvec+c
 
     # find lagrange multipliers
-    #muk = transpose(A) * Qⁱ * (A * c - q)
-    dAcdf = spzeros(length(ufvec),length(c))
-    dA1df = spzeros(length(ufvec),length(c))
+    muk = transpose(A) * Qⁱ * (A * csum - q)
+    dAcdf = spzeros(length(ufvec))#,length(c))
+    dA1df = spzeros(length(ufvec))#,length(c))
+    A1 = (A * onesvec)
+    Acq = Qⁱ*(A * csum - q)
     for ii in eachindex(ufvec)
-          dAcdf[ii,non_zero_indices[ii, 1]] = c[non_zero_indices[ii, 2]]
-          dA1df[ii,non_zero_indices[ii, 1]] = 1
+          #dAcdf[ii,non_zero_indices[ii, 1]] = c[non_zero_indices[ii, 2]]
+          #dA1df[ii,non_zero_indices[ii, 1]] = 1
+	  dAcdf[ii] = sum(csum[non_zero_indices[ii, 2]].*Acq[non_zero_indices[ii, 1]])
+	  dA1df[ii] = sum(A1[non_zero_indices[ii, 1]])
     end
-    dAdf_terms = dA1df * (A * onesvec) + dAcdf * Qⁱ * (A * c - q)# + dA1df * (A * onesvec)
+    dAdf_terms = dA1df + dAcdf #* Qⁱ * (A * c - q)# + dA1df * (A * onesvec)
 
     if guvec != nothing
         tmp = guvec
         for (ii,vv) in enumerate(tmp)
             if ii <= ulength
-               guvec[ii] = 2 * 10^3 * uvec[ii]#-(2 * transpose(muk) * Wⁱ)[ii]
+               guvec[ii] = 20 * uvec[ii]+2 * muk[ii]
             else
                guvec[ii]= 2 * dAdf_terms[ii-ulength]# + 4 * 10 * convec[ii].^3
             end
@@ -1747,9 +1755,9 @@ function costfunction_gridded_model!(J,guvec,convec::Vector{T},non_zero_indices,
     end
     
     if J !=nothing
-        return 10^3 * uvec ⋅ uvec + transpose(A * c - q) * Qⁱ * (A*c - q)+
+        return 10 * uvec ⋅ uvec + transpose(A * csum - q) * Qⁱ * (A*csum - q)+
 	#2 * transpose(muk)*( Wⁱ * uvec+c-y)+ sum(ufvec)#sum(A*onesvec)
-               transpose(A * onesvec) * (A* onesvec)# + 10 * ufvec ⋅ ufvec.^3
+                transpose(A * onesvec) * (A* onesvec)# + 10 * ufvec ⋅ ufvec.^3
     end
 end