My simple version using Image:

size = 300;
r = ListConvolve[DiskMatrix[#], 
     RandomInteger[BernoulliDistribution[0.001], {5 size, size}], {1, 1}] & /@ {1.5, 2, 3};
Dynamic[Image[(r[[#]] = RotateRight[r[[#]], #]) & /@ {1, 2, 3}; Total[r[[All, ;; size]]]]]

Update

A slightly prettier version, same basic idea but now with flakes.

flake := Module[{arm},
   arm = Accumulate[{{0, 0.1}}~Join~RandomReal[{-1, 1}, {5, 2}]];
   arm = arm.Transpose@RotationMatrix[{arm[[-1]], {0, 1}}];
   arm = arm~Join~Rest@Reverse[arm.{{-1, 0}, {0, 1}}];
   Polygon[Flatten[arm.RotationMatrix[# \[Pi]/3] & /@ Range[6], 1]]];

snowfield[flakesize_, size_, num_] := 
  Module[{x = 100/flakesize}, 
   ImageData@
    Image[Graphics[{White, 
       Table[Translate[
         Rotate[flake, RandomReal[{0, \[Pi]/6}]], {RandomReal[{0, x}],
           RandomReal[{0, 5 x}]}], {num}]}, Background -> Black, 
      PlotRange -> {{0, x}, {0, 5 x}}], ImageSize -> {size, 5 size}]];

size = 300;

r = snowfield @@@ {{1, size, 500}, {1.5, size, 250}, {2, size, 50}};
Dynamic[Image[(r[[#]] = RotateRight[r[[#]], #]) & /@ {1, 2, 3}; 
  Total[r[[All, ;; size]]]]]

I happened to create some snowflakes and snow fall a couple weeks back, and its nice to have some place to share with others!

First, we create some algorithmically generated snowflakes with some randomness using a kind of iterated function system based off the 6-pointed "star" shown below.

H = Table[{Cos[n*Pi/3], Sin[n*Pi/3]}, {n, 0, 5, 1}];

Graphics@Table[Line[{{0, 0}, H[[i]]}],{i, 1, 6, 1}]

To construct the snowflake we just make copies of this star which are translated by some random amount and rotated (for 6-fold symmetry), and scaled in size.

We can seed some random real numbers to generate unique snowflakes

rr[n_] := (SeedRandom[n]; RandomReal[])

SnowFlake[Q_, x_, y_, spin_, size_, k_, h_, opacity_, scale_, N_, PR_, IS_] :=
Graphics[{
 Rotate[
 Translate[
 Scale[
  Table[
   Table[
    Rotate[
    Translate[
    Scale[
       Table[
        {AbsoluteThickness[k*h^(n - 1)], Opacity[opacity], White,
         Line[
          {{0, 0}, H[[i]]}]},
        {i, 1, 6, 1}],   
    scale^(n - 1)],
    {If[n == 1, 0, rr[Q*n]], 0}],
    If[n == 1, 0, (j + rr[Q*n])*Pi/3], {0, 0}],
   {j, 0, 5, 1}],
  {n, 1, N, 1}],
size],
{x, y}],
spin, {x, y}]
},
PlotRange -> PR, ImageSize -> IS, Background -> Black]

There are various parameters one could adjust here including some like x, y position and spin to rotate a snowflake which will be used later to animate snowfall. Here is a sample of several snowflakes after fixing a particular choice of parameters and ranging Q to seed different random reals:

GraphicsGrid[
Table[
SnowFlake[Q*W, 0, 0, (-1)^(Round[rr[4 Q*W]]) (0 + rr[2 Q*W]) Pi/3, 1, 1,.85, .8, .5 + .2 rr[3 Q*W], 15, 2, 100],
{Q, 7, 13, 1}, {W, 14, 18, 1}], 
Background -> Black, ImageSize -> {500, 700}, AspectRatio -> 7/5]

To generate snowfall, we parameter the x,y position, and spin of each snowflake with time t. I used the following code. Notice that Snowflake[] was used twice and combined with Show[]. The reason this was done was to create a better simulation of depth in the snowfall. The only difference between the two is the size the snowflakes are allowed to randomly range through. Here, Q1 and Q2 can be varied to generate different snowflakes, and F1 and F2 control the number of snowflakes.

SF[Q1_, Q2_, F1_, F2_, t_] :=  Show[
Table[SnowFlake[
    Q1*f,
    4rr[2 Q1*f] - 2,
    1.25 - 2.5 (Mod[t + rr[3 Q1*f], 1]),
    (-1)^(Round[rr[4 Q1*f]])*t*Ceiling[6 rr[4 Q1*f]]*Pi/3,
    .04 + .03 rr[5 Q1*f],
    .8, 1, .3, .75, 6, {{-2, 2}, {-6/5, 6/5}}, {500, 300}],
   {f, 1, F1, 1}],
  Table[
   SnowFlake[
    Q2*f, 4 rr[2 Q2*f] - 2,
    1.25 - 2.5 (Mod[t + rr[3 Q2*f], 1]),
    (-1)^(Round[rr[4 Q2*f]])*t*Ceiling[6 rr[4 Q2*f]]*Pi/3,
    .07 + .12 rr[5 Q2*f],
    .8, 1, .3, .75, 6, {{-2, 2}, {-6/5, 6/5}}, {500, 300}],
   {f, 1, F2, 1}]
 ]

Export[
  "snowfall.gif",
  Table[
   SF[22, 40, 30, 15, t],
   {t, 0, 39/40, 1/40}]
  ]

This animated gif loops seamlessly with 40 frames, but can of course be made smaller if desired. In my original project, I included some text (code not shown):

This approach uses a voxel-based snowflake that actually looks quite realistic because of the way Image3D displays data. The set-up is similar to that of Simon, but everything is in 3D. So the calculations will be correspondingly more time-consuming.

flakyData = 
  ImageData[
   Dilation[
    RandomImage[BernoulliDistribution[0.0001], {150, 150, 150}], 
    DiskMatrix[{1, 1, 1}]]];

Image3D[flakyData]

With this, we now create an animation by rotating the data array downward:

AbsoluteTiming[
 snowAnim = 
   Table[Rasterize[
     Show[Image3D[flakyData = RotateRight[flakyData, 3], 
       ColorFunction -> "BlackWhiteOpacity"], 
      Background -> Darker[Gray], ViewPoint -> Right, ViewAngle -> .6,
       Boxed -> False], "Image"],
    {i, 1, 50}];]

(* ==> {14.178811, Null} *)

So you have to wait several seconds for the frames to be generated. Here is the result:

Export["snowFall.gif", snowAnim, "DisplayDurations" -> .07, 
  "AnimationRepetitions" -> Infinity];

By choosing the ViewAngle narrower than usual, the view of the 3D box in which the flakes are falling doesn't include the empty space on the sides. But you could change that to make the appearance and disappearance of the snow flakes more noticeable. That would correspond to looking out into the night with a flashlight where you see the flakes only when they enter the beam...

If you're really patient, you can add some sideways wind as Rojo was suggesting in his 3D solution, but here I would realize that with additional matrix rotations:

data = RandomInteger[BernoulliDistribution[0.0001], {150, 150, 150}];

AbsoluteTiming[
 snowAnim = 
   Table[Rasterize[
     Show[Dilation[
       Image3D[data = 
         MapIndexed[RotateRight[#, Floor[#2/30]] &, 
          RotateRight[data, 3]], 
        ColorFunction -> "BlackWhiteOpacity"], DiskMatrix[{1, 1, 1}]],
       Background -> Darker[Gray], ViewPoint -> Right, 
      ViewAngle -> .6, Boxed -> False], "Image"],
    {i, 1, 150}];]

(* ==> {142.734662, Null} *)

This is a snow storm. The reason why this takes so much longer to create is that I can't take the Dilate operation out of the Table because the RotateRight in different directions by different amounts will distort the flakes (unless they are single voxels). Maybe there's another way to speed this up, I just wanted to show what the Image3D approach looks like.

Edit: making a GIF into a dynamic animated notebook display

Since the question asked about GIF creation and Dynamic display, this is an opportunity to do both combined, with the GIF approach as described above as the first step.

The second step is now to re-import the GIF I just exported. For this, I have a general function that looks like this:

makeAnimation[list_, delayList_: {.03}] := 
 DynamicModule[{l = Length[list], delays = Abs@Flatten[{delayList}], 
   times, totalTime, delta = .03, frames},
  times = Round[.5 + PadRight[delays, l, delays]/delta];
  frames = 
   Flatten@Table[Table[list[[i]], {times[[i]]}], {i, Length[times]}];
  totalTime = Length[frames];
  EventHandler[
   Dynamic[frames[[Clock[{1, totalTime, 1}, totalTime delta]]], 
    TrackedSymbols -> {}], {"MouseUp", 2} :> Null
   ]
  ]

importGIF[fileName_] := Module[{lst, seq, pos, durations, frameList},
  lst = BinaryReadList[fileName];
  seq = FromDigits[#, 16] & /@ StringSplit["21 f9 04"];
  pos = Position[Partition[lst, 3, 1], seq];
  durations = (Extract[lst, pos + 4] + 256 Extract[lst, pos + 5])/
    100;
  frameList = Import[fileName];
  makeAnimation[frameList, durations]
  ]

importGIF["snowFall.gif"]

(* Same output as the first animation above, but now playing in notebook. *)

The output cell can be copied anywhere you like, even into a new notebook. To reiterate: you only need to copy the output by selecting its cell bracket, without having to regenerate the movie by re-evaluating the input. So the function importGIF is the closest one can get to having the original GIF play in the Mathematica notebook the way it plays in a web browser.

An important ingredient in this function is actually not really used in this context: the ability to have a varying frame rate. This utilizes a solution by Heike to the question When importing GIF animation, how to find the correct list of “DisplayDurations”? ^[1].

My attempt adresses the second option. These are simple falling circles with random shift and color depending on the radius of the "snowflake".

snow[n_, {w_, h_}] := 
  DynamicModule[{flakes, cut, offset, drawflake, rmax, mgray = 0.85},
   rmax = w/100;
   drawflake = {GrayLevel[mgray + (1 - mgray) #2], Disk[#1, #2 rmax]} &;
   offset = cut @@ (#1 + 14 {#2 RandomReal[{-0.5, 0.5}], -#2}) &;
   cut = {#1~Mod~w, #2~Mod~h} &;
   flakes = Array[{{Random[] w, Random[] h}, Random[]} &, n]~SortBy~Last;
   Dynamic[
    Graphics[
     drawflake @@@ (flakes = {offset[#1, #2], #2} & @@@ flakes),
     Background -> LightGray,
     PlotRange -> {{0, w}, {0, h}},
     ImageSize -> {w, h},
     AspectRatio -> Automatic
     ]
    ]];

The usage is

snow[200, {800, 400}]

Here 200 is the number of snowflakes and {800,400} is the dimensions of the graphics. The screenshot:

And finally a "full screen" version:

{{{z, fsw}, {z, fsh}}} = FullScreenArea /. (ScreenInformation /. 
    Options[$FrontEnd, ScreenInformation]);
CreateDocument[snow[200, {fsw, fsh}], WindowSize -> Full];

To animate snowflakes, we first need to create them. Here's a simple way to create snowflake-like objects:

makeLine = {Re@#, Im@#} &;
rot = Exp[I*Pi/3];
list = {0, 1, rot, rot^2, rot, 1, 1/rot, 1/rot^2, 1/rot^3, 1/rot^2, 
   1/rot, 1};
flake = Flatten@NestList[1/rot*# &, list, 7];
flakeprim = 
  Scale[Translate[Line[#], -Mean[#]], 1] &@
   Accumulate[makeLine /@ flake];

This uses the fact that multiplication by $\exp(i\phi)$ is a rotation by $\phi$ in the complex plane to successively rotate a straight line to create one arm of the flake, then rotating the whole thing a few times to create the whole flaks.

Here's how it looks:

Graphics[
    flakeprim,
    Axes -> True
]

Some might say that this looks more like a spaceship than a snowflake. Maybe it does.

This may be animated using a minimal modification of faleichik's code (minimal means that I left everything alone except what I had to change, so bits are redundant):

snow[n_, {w_, h_}] := 
 DynamicModule[{flakes, cut, offset, drawflake, rmax, mgray = 0.85}, 
  rmax = w/100;
  drawflake = Scale[Translate[flakeprim, #1], #2] &;
  offset = cut @@ (#1 + 14 {#2 RandomReal[{-0.5, 0.5}], -#2}) &;
  cut = {#1~Mod~w, #2~Mod~h} &;
  flakes = 
   Array[{{Random[] w, Random[] h}, Random[]} &, n]~SortBy~Last;
  Dynamic[
   Graphics[
    drawflake @@@ (flakes = {offset[#1, #2], #2} & @@@ flakes), 
    Background -> LightGray, PlotRange -> {{0, w}, {0, h}}, 
    ImageSize -> {w, h}, AspectRatio -> Automatic]]]

so eg

snow[30, {500, 500}]

This could be a start for someone who wants to play around in v9

numFlakes = 400;
flake = Sphere[{0, 0, 0}, 0.3];
wind = {2, 1};
procs = Table[
   WienerProcess[\[Mu], 1], {\[Mu], {Sequence @@ wind, -3}}];
sampleTDs = 
  Table[RandomFunction[proc, {0, 20, 0.05}, numFlakes], {proc, procs}];
funs = With[{offset = {0, 0, 20} + {20, 20, 40} RandomReal[{-1, 1}, 3]}, 
     Function[time, Through@#[time] + offset]] & /@ 
   Transpose@Through@sampleTDs["PathFunction", All];

and then

Graphics3D[{Translate[flake, 
   Dynamic[With[{time = Clock[{0, 10, 0.05}]}, #@time & /@ funs]]]},  
 Boxed -> False, Lighting -> "Neutral", AxesLabel -> {"x", "y", "z"}, 
 Background -> LightGray, 
 PlotRange -> {{-10, 10}, {-10, 10}, {0, 20}}]

I created some (semi-random) snow for a Christmas related notebook I recently wrote. Thought it was worth sharing:

FallDown[Disk[{x_, y_}, r_]] := 
 Disk[{x + RandomReal[{-0.03, 0.05}], y - RandomReal[{0.05, 0.25}]}, r]
Snow[] := Module[{
   snowflakes = Disk[#, RandomReal[{0.03, 0.1}]] & /@
     Flatten[
      Table[{i + RandomReal[{-0.3, 0.3}], 
        j + RandomReal[{-0.3, 0.3}]}, {i, 0, 12}, {j, -5, 20}], 1],
   list = {}},
  Do[
   snowflakes = FallDown /@ snowflakes;
   AppendTo[list, Graphics[{White, snowflakes},
     PlotRange -> {{0, 12}, {0, 12}},
     Background -> Black]],
   {52}];
  list]

The final output of the notebook also adds some (again semi-randomly placed) trees and hills, along with Santa himself :-)

You may be wondering where all the snow seen in the brilliant and amazing answers above (assuming you've sorted answers by votes) ends up... The answer is, down here, on the floor. Sadly, I can't upload this riveting movie in its full unabridged splendour, although at least my cat enjoys watching the full 3-hour version.

This is how it starts:

and here is the epic finale:

flakes1 = 
  Table[Text[
    Style[FromCharacterCode[16^^2744], RandomInteger[{24, 30}], 
     White], {RandomReal[], RandomReal[]}], {x, 3000}]; 
flakes2 = 
  Table[Text[
    Style[FromCharacterCode[16^^2746], RandomInteger[{24, 30}], 
     White], {RandomReal[], RandomReal[]}], {x, 3000}] ;
flakes = Riffle[flakes1, flakes2];
t = Table[
   Graphics[
    {
     Opacity[0.8],
     flakes[[1 ;; x]]
     }, 
     Background -> Gray, 
     PlotRange -> {{0, 1}, {0, 1}}],
   {x, 1, Length[flakes]}];

Part of the answer lies in generating snowflakes.

If one is truly interested in a random-generated snowflake, I suggest using varying initial conditions and rules for a 2D Cellular Automaton on a hexagonal grid to generate them. A demonstration for this already exists. ^[1]