def METHOD_NAME(data, label_attr=None):
    ana_data = {}
    xlim = {}
    # plot instantaneous profiles
    varlist = ['salt_3d', 'buoy_freq_3d', 'shear_freq_3d', 'tke_3d', 'eps_3d', 'len_3d', 'eddy_diff_3d']
    log_variables = []
    nplots = len(varlist)
    fig, axlist = plt.subplots(nrows=1, ncols=nplots, sharey=True,
                               figsize=(nplots*2.3, 6))
    # plot analytical solutions
    for v, ax in zip(varlist, axlist):
        if v in ana_data:
            zz, uu = ana_data[v]
            ax.plot(uu, zz, 'r', lw=1.7, linestyle='dashed',
                    label='analytical', zorder=10)
        ax.grid(True)
        ax.set_title(field_metadata[v]['shortname'].replace(' ', '\n'))
        ax.set_xlabel('[{:}]'.format(field_metadata[v]['unit']), horizontalalignment='right')
        loc = matplotlib.ticker.MaxNLocator(nbins=3, prune='upper')
        fmt = matplotlib.ticker.ScalarFormatter(useOffset=None, useMathText=None, useLocale=None)
        fmt.set_powerlimits((-2, 3))
        ax.xaxis.set_major_locator(loc)
        ax.xaxis.set_major_formatter(fmt)
        if v in xlim:
            ax.set_xlim(xlim[v])
    # plot_models
    for entry in data:
        for v, ax in zip(varlist, axlist):
            solver_obj = entry['solver']
            layers = entry['nz']
            z, xyz = construct_plot_coordinates(layers)
            func = solver_obj.fields[v]
            arr = numpy.array(func.at(tuple(xyz)))
            print('field: {:} min {:} max {:}'.format(v, arr.min(), arr.max()))
            if len(arr.shape) == 2:
                # take first component of vectors
                arr = arr[:, 0]
            label = model_label(entry, label_attr=label_attr)
            if v in log_variables:
                ax.semilogx(arr, z, lw=1.5, ls='solid', label=label, alpha=0.7)
            else:
                ax.plot(arr, z, lw=1.5, ls='solid', label=label, alpha=0.7)
    axlist[0].set_ylabel('z [m]')
    axlist[0].set_ylim([-depth*1.005, 0])
    # add legend
    plt.legend(loc='upper left', bbox_to_anchor=(1.02, 1.0),
               prop={'size': 'small'})
    closure_str = '-'.join(sorted(unique([closure_abbr[e['gls_closure']].replace('-', '') for e in data])))
    stab_str = '-'.join(sorted(unique([stab_abbr[e['stability_func']] for e in data])))
    nz_str = '-'.join(map(str, sorted(unique([e['nz'] for e in data]))))
    imgfile = 'kato_profiles_{:}_{:}_nz{:}.png'.format(
        closure_str, stab_str, nz_str)
    print('Saving figure {:}'.format(imgfile))
    plt.savefig(imgfile, bbox_inches='tight', dpi=200.)
    plt.close(fig)