For the kriging step, you should probably use Empirical Bayesian Kriging (EBK). It is our most accurate interpolation method, and it is easy to automate because it is available as a geoprocessing tool (it can also be accessed from the Geostatistical Wizard). You'll want to output a raster and set both the Extent and Mask environments to your boundary shape file. You can also use the Output Coordinate System environment in EBK to define the spatial reference of the output raster.
I'm not sure how to do the rest, but hopefully someone else can help with those parts.