raster - using Python to calculate NDVI, output as a Geotiff

I need a script that can calculate NDVI from two separate input Geotiff files, then output the results as a Geotiff. I wrote following script, most of which is borrowed from various on-line sources (I've never had any training with python and only minimal experience). I know the calculation part works, because if I print out the ndvi array, I get reasonable values. If I calculate statistics on the array, I also get reasonable results. The script does create an output file called ndvi.tif, but all the values in the file are 0 ( and not 0.0).

I need to convert this into a batch script that can do this calculate on a large number of input files, but right now I can't figure out why it won't output calculated values to a file. Any advice?

 +++++++++++++++++++++++++++++++++++++++++++++

 from osgeo import gdal
 import numpy as np

 from numpy import *
 g = gdal.Open("fhmod1.tif")
 red = g.ReadAsArray()
 g = gdal.Open("fhmod2.tif")
 nir = g.ReadAsArray()
 red = array(red, dtype = float)
 nir = array(nir, dtype = float)
 check = np.logical_and ( red > 1, nir > 1 )
 ndvi = np.where ( check,  (nir - red ) / ( nir + red ), -999 ) 
 geo = g.GetGeoTransform()  

 proj = g.GetProjection()   
 shape = red.shape        
 driver = gdal.GetDriverByName("GTiff")
 dst_ds = driver.Create( "ndvi.tif", shape[1], shape[0], 1, gdal.GDT_Float32)
 dst_ds.SetGeoTransform( geo )
 dst_ds.SetProjection( proj ) 
 dst_ds.GetRasterBand(1).WriteArray(ndvi)