Geiger's method of locating local earthquakes has been extended to include the effect of P velocity variation along the ray paths in three dimensions. The crustal structure was modeled by rectangular blocks, and a parameter was assigned to each block describing the perturbation of P wave slowness in the block. On the basis of an initial model, a set of linear equations for the observed first P arrival times was formulated in terms of the source and medium parameters. The source parameters for all the earthquakes in the data set and medium parameters for all the blocks penetrated by the seismic rays were then determined simultaneously by the damped least squares method. A computer program has been written for the case of a homogeneous initial medium model with constant P velocity. It has been tested with artificial data and applied to a set of a local earthquakes recorded by a dense seismic array in Bear Valley, California. The resultant velocity distribution in the top 5 km shows a narrow low-velocity zone of about 5km/s in the San Andreas fault zone sandwiched between high-velocity regions of about 6km/s.