model_dict.py

# Start by storing model names in an array. This keeps order, but prolly not important
models = [ 'casa', 'lpj', 'jules', 'middle_aged_glamour_model']

# The year list that you are going to loop over
years = np.arange ( 1990, 2011 ) 

# Model variables for each model
# You need to keep each list in the order they appear in the output. That's fairly important.
model_var = {  'casa':[ 'doy','NEE', 'NPP'], 'lpj':['year', 'doy','NEE', "GPP"], \
                       'doy',middle_aged_glamour_model']=[ 'doy','you','dont','wanna_know'] }
# Initialise the output "object". It's an empty dictionary
# if you need multiple treatments, define two model outputs and add an extra external loop.
model_outputs = {}
for model in models: # loop over models
  model_outputs [ model ] = {} 
  # The container will be a dictionary with keys given by model names. Nice for plot names and stuff
  for year in years:   # loop over years
       # Each year is stored as its own key. If years are integers, else convert into integers
       model_outputs [ model ][ year ] = {} # empty dictionary again...
       retval = run_model ( model, year ) # Model run, retval is a tuple or an array or whatever
       i = 2 # retval [0] is year, retval[1] is doys, retval[2]... retval[N] is data arrays
       # Loop over model variables (remember they are in order?) and add dictionary entries
       for vari in model_var [ model ]:
          # now the entry in the dictionary is an array (usually of length 365) with the parameters
          # of interest. I am assuming that retval looks like this
          # retval = ( year, doy_arr, param1_arr, param2_arr, ... )
          # year is an integer
          # doy_arr is a numpy array of eg. int16
          # param1_arrX are arrays of the same shape as doy_arr storing the daily values of the parameters
          model_outputs[ model][year][vari] = retval[i]
          i = i + 1
       model_outputs[ model][year]['doy'] = retval[1]

# So you can look at the NPP from LPJ in 2004 as model_outputs['lpj'][2004]['NPP'], a numpy array

model_var - i.e. the first bit that contains the headers varies with models so needs to go within the loop

model_var is defined per model (it is a dictionary, keys being the model name, items stored under that key are a list of model outputs).

OK almost have it though I am getting junk out...

def grab_all_model_data(self):

    models = ['CASA', 'CLMX', 'EALC', 'EDXX', 'GDAY', 'LPJG', 'ISAM', 'OCNX', 'SDVM', 'TECO'] 
    years = np.arange(1990, 2011)

    model_var = {}
    model_outputs = {}
    for model in self.model_dict:
        where = self.model_dir + model

        model_outputs[model] = {}
        # read header info and then open the data
        fname = os.path.join(where, 'Day_' + model + self.site + 'AMB' + '.asc.hdr')
        (header, numrows, numcols) = self.read_header(fname)
        model_var[model] = header


        #model_outputs[model] = treatment
        fname = os.path.join(where, 'Day_' + model + self.site + 'AMB' + '.asc')
        data = gd.fillUpArray(fname).reshape(numrows, numcols)
        i = 0
        for var in model_var[model]:
            model_outputs[model][var] = data[i]
            i += 1


    print(model_outputs['CASA']['YEAR'])

However whilst I would like to have that year loop like you do - not all the models start and end in the same point. So ideally I would like to fill with undefined to make sure it doesn't crash later

maaaaaaaaaaaaaannnn relax!! I have set it to be a dictionary and year is part of the key, so you can check whether a model has output for a given year, and skip chucking it into the dictionary if it doesn't. This is what you get from learning C and awk and gnuplot... now you have to un-learn it for the 21st century :-)