Adding time dimension and corresponding new variable to a netcdf file
There are two issues (and one optimization) with/for the code above:


Array indexing starts with 0 (not 1) in Python. Therefore, it should be if data[0,m] == countries[i,j]: because you are interested in the first row;
m goes from 0 to 1 (for m in xrange(0,1):) but it should go from 0 to number_of_countries so that all columns of the energy are iterated. 
The order of the n and m for loops seems to be quite ineffective. Better would be: for m ...: if data[1,m] == ...: for n in ...: nc.variables...=.... Or replace latter for loop by nc.variables["energy"][:,i,j]=data[:,m]


This could work (please improve => community Wiki) but I cannot test it due to missing input data:

y = nc.createDimension("y", 950)  #y
x = nc.createDimension("x", 1000) #x
time = nc.createDimension("time", 87648) #time

latitude=nc.createVariable("y", "f8", ("y", ))
longitude=nc.createVariable("x", "f8", ("x", ))
time = nc.createVariable("time","f8",("time",))
countries = nc.createVariable("countries", "f8", ("y", "x"));
energy = nc.createVariable("energy", "f8", ("time", "y", "x"), 
                           fill_value=-9999, chunksizes=(1, 950, 1000));

time.standard_name='time'
time.units ='hours since 2006-01-01 00:00:00.0'
time.calendar='proleptic_gregorian'

lats =  np.arange(5497500,747500,-5000)    #y  
lons =  np.arange(2502500,7502500,5000)    #x 

# Fill variables
latitude[:] = lats
longitude[:] = lons

countries [:]=value

for i in xrange(0,950):  # grid rows
    for j in xrange(0,1000):    # grid columns
        for m in xrange(0,int(max(countries))):    # countries
            if data[0,m] == countries[i,j]:   #correspondance
                nc.variables["energy"][:,i,j]=data[1:87648,m]