I have a particle where I know:

1. position (either in euclidean or spherical coordinates)
2. direction (360°/2 pi radians)
3. speed (linear distance).

How do one calculate next position for the particle?

Weird. I just answered a question about quaternions that is essentially the same question.

If your time step is small, move the particle in R^3, forgetting about the sphere for a moment, then project back to the sphere. If that approximation is not good enough, there is a formula you can use that involves computing cosine and sine. I'll figure it out and post it, if you really need it.

I don't know how to calculate forward in R^3 either. One solution I thought of is to rotate particle direction around sphere normal at particle position. Sound simple but I find matrix rotations hard. I don't understand how to write a "vector3.rotateAround(vector axis, float angle)"-function.

My velocities are probably small but I don't know yet, I trying to simulate fluids and they have a tendency to blow up if not correct.

Right now I would settle for any algorithm that would help me progress.

If (x,y,z) is a point on the surface of the unit sphere centered at the origin, the East vector is something like (z, 0, -x) and the North vector is something like (-xy/S, S, -yz/S), with S:=sqrt(1-y^2). You can then express your direction as North * cos(angle) + East * sin(angle).

I kind of rushed through the Math, so I am not 100% certain that the formulas are correct, but you can try to reproduce them yourself.

I don't understand how to write a "vector3.rotateAround(vector axis, float angle)"-function.

It's very easy to do with quaternions. You just need to know how to construct the quaternion that represents the rotation you want, which is cos(angle/2) + sin(angle/2)*x*i + sin(angle/2)*y*j + sin(angle/2)*z*k, and how to rotate a point by the rotation represented by a unit-length quaternion, which is q' = q * v * conj(q).

If (x,y,z) is a point on the surface of the unit sphere centered at the origin, the East vector is something like (z, 0, -x) and the North vector is something like (-xy/S, S, -yz/S), with S:=sqrt(1-y^2). You can then express your direction as North * cos(angle) + East * sin(angle).

I kind of rushed through the Math, so I am not 100% certain that the formulas are correct, but you can try to reproduce them yourself.

Really, thanks a lot!

I tried it and it works.