import pyaerocom as pya
def read_modeldata_pyaerocom(MODEL_ID, VAR, REG):
    pya.const.BASEDIR = '/home/notebook/shared-ns1000k/inputs/pyaerocom-testdata/'
    reader = pya.io.ReadGridded(MODEL_ID)
    VAR = reader.read_var(VAR, start=2010, ts_type='monthly')
    if REG=='ARCTIC':
        VAR = VAR.crop(lat_range=(70,90))
    df = VAR.crop(lat_range=(70,90)).to_xarray().to_dataframe()
    return df

MODEL_ID = 'ECHAM6.3-HAM2.3-met2010_AP3-CTRL'
VAR = [['od550bc','od550so4','od550du','od550ss'],['loadbc','loadso4','loaddu','loadss']]
REG = 'ARCTIC'
VAR_MEAN = [[],[]]
for i in range(len(VAR)):
    for j in range(len(VAR[i])):
        df = read_modeldata_pyaerocom(MODEL_ID, VAR[i][j], REG)
        vm = df[VAR[i][j]].mean()
        VAR_MEAN[i].append(vm)

MEC_MEAN = []
for i in range(len(VAR[0])):
    MEC_MEAN.append(VAR_MEAN[0][i]/VAR_MEAN[1][i])

MEC_MEAN