Is it possible to rotate an object around its own axis and not around the base coordinate's axis?

First thing you should know is that in OpenGL, transformation matrices are multiplied from right. What does it mean? It means that the last transformation you write gets applied to the object first.

So let's look at your code:

gl.glScalef(0.8f, 0.8f, 0.8f);
gl.glTranslatef(0.0f, 0.0f, -z);
gl.glRotatef(xrot, 1.0f, 0.0f, 0.0f);   //X
gl.glRotatef(yrot, 0.0f, 1.0f, 0.0f);   //Y
gl.glRotatef(zrot, 0.0f, 0.0f, 1.0f);   //Z
gl.glTranslatef(0.0f, 0.0f, z);

square.draw(gl);    

This means that, first, the object is moved to (0.0f, 0.0f, z). Then it is rotated around Z, then around Y, then around X, then moved by (0.0f, 0.0f, -z) and finally scaled.

You got the scaling right. You put it first, so it gets applied last. You also got

gl.glTranslatef(0.0f, 0.0f, -z);

in the right place, because you first want to rotate the object then move it. Note that, when you rotate an object, it ALWAYS rotates around the base coordinate, that is (0, 0, 0). If you want to rotate the object around its own axes, the object itself should be in (0, 0, 0).

So, right before you write

square.draw(gl);

you should have the rotations. The way your code is right now, you move the object far (by writing

gl.glTranslatef(0.0f, 0.0f, z);

before square.draw(gl);) and THEN rotate which messes things up. Removing that line gets you much closer to what you need. So, your code will look like this:

gl.glScalef(0.8f, 0.8f, 0.8f);
gl.glTranslatef(0.0f, 0.0f, -z);
gl.glRotatef(xrot, 1.0f, 0.0f, 0.0f);   //X
gl.glRotatef(yrot, 0.0f, 1.0f, 0.0f);   //Y
gl.glRotatef(zrot, 0.0f, 0.0f, 1.0f);   //Z

square.draw(gl);    

Now the square should rotate in place.

Note: After you run this, you will see that the rotation of the square would be rather awkward. For example, if you rotate around z by 90 degrees, then rotating around x would look like rotating around y because of the previous rotation. For now, this may be ok for you, but if you want to it to look really good, you should do it like this:

Imagine, you are not rotating the object, but rotating a camera around the object, looking at the object. By changing xrot, yrot and zrot, you are moving the camera on a sphere around the object. Then, once finding out the location of the camera, you could either do the math and get the correct parameters to call glRotatef and glTranslatef or, use gluLookAt.

This requires some understanding of math and 3d imagination. So if you don't get it right the first day, don't get frustrated.

Edit: This is the idea of how to rotate along rotated object coordinates;

First, let's say you do the rotation around z. Therefore you have

gl.glRotatef(zrot, 0.0f, 0.0f, 1.0f);   //Z

Now, the global Y unit vector is obviously (0, 1, 0), but the object has rotated and thus its Y unit vector has also rotated. This vector is given by:

[cos(zrot) -sin(zrot) 0]   [0]   [-sin(zrot)]
[sin(zrot)  cos(zrot) 0] x [1] = [ cos(zrot)]
[0          0         1]   [0]   [ 0        ]

Therefore, your rotation around y, should be like this:

gl.glRotatef(yrot, -sin(zrot), cos(zrot), 0.0f);   //Y-object

You can try this so far (disable rotation around x) and see that it looks like the way you want it (I did it, and it worked).

Now for x, it gets very complicated. Why? Because, the X unit vector is not only first rotated around the z vector, but after it is rotated around the (-sin(zrot), cos(zrot), 0) vector.

So now the X unit vector in the object's cooridnate is

                   [cos(zrot) -sin(zrot) 0]   [1]                      [cos(zrot)]
Rot_around_new_y * [sin(zrot)  cos(zrot) 0] x [0] = Rot_around_new_y * [sin(zrot)]
                   [0          0         1]   [0]                      [0        ]

Let's call this vector (u_x, u_y, u_z). Then your final rotation (the one around X), would be like this:

gl.glRotatef(xrot, u_x, u_y, u_z);   //X-object

So! How to find the matrix Rot_around_new_y? See here about rotation around arbitrary axis. Go to section 6.2, the first matrix, get the 3*3 sub matrix rotation (that is ignore the rightmost column which is related to translation) and put (-sin(zrot), cos(zrot), 0) as the (u, v, w) axis and theta as yrot.

I won't do the math here because it requires a lot of effort and eventually I'm going to make a mistake somewhere around there anyway. However, if you are very careful and ready to double check them a couple of times, you could write it down and do the matrix multiplications.

Additional note: one way to calculate Rot_around_new_y could also be using Quaternions. A quaternion is defined as a 4d vector [xs, ys, zs, c], which corresponds to rotation around [x, y, z] by an angle whose sin is s and whose cos is c.

This [x, y, z] is our "new Y", i.e. [-sin(zrot), cos(zrot), 0]. The angle is yrot. The quaternion for rotation around Y is thus given as:

q_Y = [-sin(zrot)*sin(yrot), cos(zrot)*sin(yrot), 0, cos(yrot)]

Finally, if you have a quaternion [a, b, c, d], the corresponding rotation matrix is given as:

[1 - 2b^2 - 2c^2        2ab + 2cd           2ac - 2bd   ]
[   2ab - 2cd        1 - 2a^2 - 2c^2        2bc - 2ad   ]
[   2ac - 2bd           2bc + 2ad        1 - 2a^2 - 2b^2]

I know next-to-nothing about openGL, but I imagine translating to 0, rotating and then translating back should work...

gl.glTranslatef(-x, -y, -z);
gl.glRotatef(xrot, 1.0f, 0.0f, 0.0f);   //X
gl.glRotatef(yrot, 0.0f, 1.0f, 0.0f);   //Y
gl.glRotatef(zrot, 0.0f, 0.0f, 1.0f);   //Z
gl.glTranslatef(x, y, z);