ListContourPlot Interpolation Order

I am trying the use ListContourPlot with interpolation of the measured data to get smooth contour lines. I tried the methods in Data interpolation and ListContourPlot but they work only for some data. I tried also with Tecplot 360 EX, where it works very well with "Kriging" however in Mathematica 12.2 I couldn’t figure it out. Any suggestion will be precious.

my data

    data={{103.8, -89.7, 27.6491}, {69.4, -56.8, 38.178}, {49.9, -30.1, 
  51.1585}, {28.4, -26.7, 54.3483}, {15, -27.1, 53.5274}, {16.3, -8.4,
   54.7628}, {166.1, -138.6, 21.191}, {12.3, -137.4, 
  26.4762}, {11.9, -114.8, 22.4272}, {6.1, -2.9, 
  60.9066}, {4.3, -14.2, 56.4494}, {0, -75, 39.1788}, {31.94, -31.81, 
  59.4118}, {0, -25, 55.0237}, {0, -40, 50.2579}, {89.7, -103.8, 
  27.6491}, {56.8, -69.4, 38.178}, {30.1, -49.9, 
  51.1585}, {26.7, -28.4, 54.3483}, {27.1, -15, 53.5274}, {8.4, -16.3,
   54.7628}, {138.6, -166.1, 21.191}, {137.4, -12.3, 
  26.4762}, {114.8, -11.9, 22.4272}, {2.9, -6.1, 
  60.9066}, {14.2, -4.3, 56.4494}, {75, 0, 39.1788}, {31.81, -31.94, 
  59.4118}, {25, 0, 55.0237}, {40, 0, 50.2579}, {5.54, 0., 
  68.5401}, {5.26885, -1.71195, 68.5401}, {4.48195, -3.25633, 
  68.5401}, {3.25633, -4.48195, 68.5401}, {1.71195, -5.26885, 
  68.5401}, {3.39227*10^-16, -5.54, 68.5401}, {40.8953, -35.3553, 
  68.4351}, {40.6242, -37.0673, 68.4351}, {39.8373, -38.6117, 
  68.4351}, {38.6117, -39.8373, 68.4351}, {37.0673, -40.6242, 
  68.4351}, {35.3553, -40.8953, 68.4351}, {33.6434, -40.6242, 
  68.4351}, {32.099, -39.8373, 68.4351}, {30.8734, -38.6117, 
  68.4351}, {30.0865, -37.0673, 68.4351}, {29.8153, -35.3553, 
  68.4351}, {30.0865, -33.6434, 68.4351}, {30.8734, -32.099, 
  68.4351}, {32.099, -30.8734, 68.4351}, {33.6434, -30.0865, 
  68.4351}, {35.3553, -29.8153, 68.4351}, {37.0673, -30.0865, 
  68.4351}, {38.6117, -30.8734, 68.4351}, {39.8373, -32.099, 
  68.4351}, {40.6242, -33.6434, 68.4351}, {35.3553, -35.3553, 
  68.4351}}

and my simple code :

    ListContourPlot[data, InterpolationOrder -> 3, Mesh -> 50, 
 MeshFunctions -> {#1 &, #2 &}, MeshStyle -> {Dotted, Dotted}, 
 Contours -> 20, ContourLabels -> All, 
 ColorFunction -> "TemperatureMap", Method -> "Kriging", 
 PlotRange -> {10, 70}]

and result

data2=Drop[data,{37}]/.{x_,y_,z_}:>{{x,y},z};
f=ResourceFunction["PolyharmonicSplineInterpolation"][data2,
   Compiled->True, InterpolationOrder->3];
Plot3D[f[x,y],{x,20,50},{y,-50,-20}]

ContourPlot[f[x,y],{x,20,50},{y,-50,-20}]