You can simply incorporate the Spatial Join tool, which is available in ArcGIS 9.3, into your model.
Alternatively, you can use python scripting in 9.3:
# Create the geoprocessor object
import arcgisscripting, sys
gp = arcgisscripting.create()
# Want to join USA cities to states and calculate the mean city population
# for each state
targetFeatures = "C:/data/USA.gdb/states"
joinFeatures = "C:/data/USA.gdb/cities"
# Output will be the target features, states, with a mean city population field (mcp)
outfc = "C:/data/USA.gdb/states_mcp"
# Create a new fieldmappings and add the two input feature classes.
fieldmappings = gp.CreateObject("FieldMappings")
# First get the POP1990 fieldmap. POP1990 is a field in the cities feature class.
# The output will have the states with the attributes of the cities. Setting the
# field's merge rule to mean will aggregate the values for all of the cities for
# each state into an average value. The field is also renamed to be more appropriate
# for the output.
fieldmap = fieldmappings.GetFieldMap(fieldmappings.FindFieldMapIndex("POP1990"))
# Get the output field's properties as a field object
field = fieldmap.OutputField
# Rename the field and pass the updated field object back into the field map
field.Name = "mean_city_pop"
field.AliasName = "mean_city_pop"
fieldmap.OutputField = field
# Set the merge rule to mean and then replace the old fieldmap in the mappings object
# with the updated one
fieldmap.MergeRule = "mean"
# Delete fields that are no longer applicable, such as city CITY_NAME and CITY_FIPS
# as only the first value will be used by default
x = fieldmappings.findfieldmapindex("CITY_NAME")
x = fieldmappings.findfieldmapindex("CITY_FIPS")
#Run the Spatial Join tool, using the defaults for the join operation and join type
gp.SpatialJoin(targetFeatures, joinFeatures, outfc, "#", "#", fieldmappings)