def interpolate_tempfilt_loop(tempfilt, zgrid, zi, output):
    """    
    Linear interpolate an Eazy "tempfilt" grid at z=zi.  
    
    `tempfilt` is [NFILT, NTEMP, NZ] integrated flux matrix
    `zgrid` is [NZ] redshift grid
    `output` is empty [NFILT, NTEMP] grid to speed up execution
    """
    sh = tempfilt.shape
    NF, NT, NZ = sh[0], sh[1], sh[2]
    #output = np.zeros((NF, NT))
    for iz in range(NZ-1):
        dz = zgrid[iz+1]-zgrid[iz]
        fint = 1 - (zi-zgrid[iz])/dz
        if (fint > 0) & (fint <= 1):
            fint2 = 1 - (zgrid[iz+1]-zi)/dz
            # print iz, zgrid[iz], fint, fint2
            for ifilt in range(NF):
                for itemp in range(NT):
                    #print ifilt, itemp
                    output[ifilt, itemp] = tempfilt[ifilt, itemp, iz]*fint + tempfilt[ifilt, itemp, iz+1]*fint2
            break              
    
    return output