# 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