Fitting a 2D elliptic paraboloid

Question

I have a set of data as follows:

data={{-0.101, 1/100000, 4.88376*10^-8}, {-0.102, 1/100000, 
  6.71978*10^-8}, {-0.1025, 1/100000, 7.75107*10^-8}, {-0.1025, 1/
  500000, 1.01625*10^-9}, {-0.1035, 1/500000, 
  4.64517*10^-10}, {-0.104, 1/500000, 1.42871*10^-9}, {-0.1055, 1/
  500000, 9.36241*10^-9}, {-0.106, 1/500000, 1.36521*10^-8}, {-0.1025,
   1/1000000, 4.25775*10^-9}, {-0.105, 1/500000, 
  5.91737*10^-9}, {-0.105, 1/1000000, 
  1.45484*10^-9}, {-0.105, -(1/1000000), 1.49644*10^-9}, {-0.105, 
  5.*10^-8, 3.71906*10^-12}, {-0.105, -5.*10^-8, 
  4.24029*10^-11}, {-0.104, -(1/1000000), 
  5.95057*10^-9}, {-0.104, -5.*10^-8, 1.80357*10^-9}, {-0.104, 0, 
  1.6742*10^-9}, {-0.104, 5.*10^-8, 1.52523*10^-9}, {-0.104, 1/
  1000000, 7.90351*10^-11}, {-0.1045, -(1/1000000), 
  3.38644*10^-9}, {-0.1045, -5.*10^-8, 5.08869*10^-10}, {-0.1045, 0, 
  4.3137*10^-10}, {-0.1045, 5.*10^-8, 3.58672*10^-10}, {-0.1045, 1/
  1000000, 3.49055*10^-10}, {-0.1051, -(1/1000000), 
  1.2098*10^-9}, {-0.1051, -5.*10^-8, 3.24277*10^-11}, {-0.1051, 0, 
  1.6652*10^-11}, {-0.1051, 5.*10^-8, 3.55996*10^-11}, {-0.1051, 1/
  1000000, 1.77581*10^-9}, {-0.1055, -(1/1000000), 
  3.9831*10^-10}, {-0.1055, -5.*10^-8, 3.89284*10^-10}, {-0.1055, 0, 
  4.13136*10^-10}, {-0.1055, 5.*10^-8, 4.94171*10^-10}, {-0.1055, 1/
  1000000, 3.38768*10^-9}, {-0.106, -(1/1000000), 
  1.32069*10^-10}, {-0.106, -5.*10^-8, 1.58681*10^-9}, {-0.106, 0, 
  1.66332*10^-9}, {-0.106, 5.*10^-8, 1.81924*10^-9}, {-0.106, 1/
  1000000, 6.15182*10^-9}}

Which was generated by simulation studies by varying the first 2 points of the data. The 3rd point in this case is the reduced chisquare when compared to a reference data with x=-0.105 and y=0. Since the data generation takes quite a bit of time usually a day or more for one single point, I can't do a normal fit on it. Therefore, I try to generate the data and do a minimized chi square fit by fitting a parabola to the data. Up till now, this approach worked as I only had 1D data (x varying over -0.104 to -0.106 or whatever value) but now with the addition of the second dimension, I am unsure how to do a reduced chisquare fit on this data. The error bars on the reduced chisquare are given by:

err={4.54114*10^-10, 5.47753*10^-10, 5.95075*10^-10, 9.22271*10^-11, 
 4.10206*10^-11, 7.68015*10^-11, 2.23605*10^-10, 2.74537*10^-10, 
 1.69164*10^-10, 1.74142*10^-10, 8.62569*10^-11, 8.59748*10^-11, 
 4.25884*10^-12, 1.86321*10^-11, 2.51042*10^-10, 1.41732*10^-10, 
 1.37909*10^-10, 1.32204*10^-10, 3.31874*10^-11, 1.85544*10^-10, 
 7.39792*10^-11, 6.80882*10^-11, 6.24258*10^-11, 5.25073*10^-11, 
 1.05285*10^-10, 2.77994*10^-11, 1.39045*10^-11, 1.96538*10^-11, 
 1.32507*10^-10, 5.69372*10^-11, 7.24441*10^-11, 6.94564*10^-11, 
 7.5394*10^-11, 1.87443*10^-10, 4.62347*10^-11, 1.40189*10^-10, 
 1.3964*10^-10, 1.45333*10^-10, 2.56496*10^-10}

I have tried fitting a fitting function:

fitModel[a_, par_, a00_, par00_, acons_, pcons_, zcons_, th_] := 
 zcons (((a - a00) Cos[th] - (par - par00) Sin[th])^2/
    acons^2 + ((a - a00) Cos[th] + (par - par00) Sin[th])^2/pcons^2)

Where a and par are known variables. However, from this I get results that I can't really trust as the fitting returns something I can't compare to the provided data. I was wondering if someone had a better solution to fitting this to find the correct fitted value (a0 and par0 in my case (which from theory I know to be -0.105 and 0 for a00 and par0 respectively)). I am aware that in my case, there will be strong correlations with all the other variables (th, offset and yscal) and that is something which can not be avoided due to the time it takes to generate the data under normal conditions.

Also I am unsure which method in the Nonlinearmodelfit is likely to result in a better result.