diff --git a/doc/userdoc.org b/doc/userdoc.org
index a88ad5da03a0c73a06cebf3e177013241cee2c8e..c92563d2c5d6fb67fab8abf1ee993087bff16f15 100644
--- a/doc/userdoc.org
+++ b/doc/userdoc.org
@@ -4106,8 +4106,8 @@ Sets verbosity mode to ~true~ or ~false~ for linear solvers.
 
 *** The ~writer~ module
 
-This module provides functionalities to numerical results to files for
-post processing.
+This module provides functionalities to write numerical results to
+files for post processing.
 
 It provides the following functions and types.
 #+NAME: get-module-info-writer
@@ -4116,7 +4116,610 @@ It provides the following functions and types.
 #+END_SRC
 #+RESULTS: get-module-info-writer
 
+**** ~writer~ provided types
+
+***** ~output~
+
+This type is used to describe output of discrete functions (of type
+~Vh~). Mainly, it associates a name to the discrete function.
+
+#+BEGIN_note
+At writing, ~pugs~ checks that the names in an ~output~ list are all
+different.
+#+END_note
+
+***** ~writer~
+
+Variables of this type manage outputs: which format is used and
+eventually writing policy. This policy sets for instance the time
+period of for time dependent post processing.
+
+**** ~writer~ provided functions
+
+***** ~name_output: Vh*string -> output~
+
+This function give a name to a discrete function.
+
+#+BEGIN_SRC pugs :exports both :results none
+  import mesh;
+  import scheme;
+  import writer;
+
+  let m:mesh, m = cartesianMesh([0], [1], 100);
+
+  let identity:R^1 -> R, x -> x[0];
+  let xh:Vh, xh = interpolate(m, P0(), identity);
+
+  let x2h:Vh, x2h = xh*xh;
+
+  write(gnuplot_writer("writer-example1"), name_output(x2h, "x_square"));
+#+END_SRC
+
+See section [[gnuplot-writers]] for details about these writers family.
+
+#+NAME: writer-example1-img
+#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/writer-example1.png")
+  reset
+  set grid
+  set border
+  unset key
+  set xtics
+  set ytics
+  set square
+  set terminal png truecolor enhanced size 640,640
+  plot '<(sed "" $PUGS_SOURCE_DIR/doc/writer-example1.gnu)' lw 2 w l
+#+END_SRC
+
+#+CAPTION: Simple illustration of the output.
+#+NAME: fig:writer-example1
+#+ATTR_LATEX: :width 0.38\textwidth
+#+ATTR_HTML: :width 300px;
+#+RESULTS: writer-example1-img
+
+
+#+BEGIN_warning
+It can be convenient to define variables or lists of ~output~. However,
+one has to pay attention of the following fact. When declaring an
+output variable, a reference *to the data* of the discrete function is
+done. It is not a reference to the language variable (or
+expression).
+#+END_warning
+
+Let us illustrate it by an important second example.
+
+#+BEGIN_SRC pugs :exports both :results none
+  import mesh;
+  import scheme;
+  import writer;
+  import math;
+
+  let m:mesh, m = cartesianMesh([0], [1], 100);
+
+  let identity:R^1 -> R, x -> x[0];
+
+  let fh:Vh, fh = interpolate(m, P0(), identity);
+  let gh:Vh, gh = fh*fh;
+
+  let pi:R, pi = acos(-1);
+
+  let output_list:(output),
+      output_list = (name_output(fh, "f"),
+                     name_output(gh, "g"),
+                     name_output(sin(pi*fh), "sin"));
+
+  let abs_cos:R^1 -> R, x -> abs(cos(pi*x[0]));
+
+  fh = interpolate(m, P0(), abs_cos); // fh now refers to another function
+
+  write(gnuplot_writer("writer-example2"), output_list);
+
+  output_list = (name_output(fh, "f"),
+                 name_output(gh, "g"));
+  write(gnuplot_writer("writer-example2-2"), output_list);
+#+END_SRC
+
+Running this code produces the gnuplot file displayed on Figure
+[[fig:writer-example2]]. One sees that ~f~ is the $\mathbb{P}_0(\mathbb{R})$
+function corresponding to the function $x \to x$ and not to the function
+$x -> |\cos(\pi x)|$. This later function is plotted on Figure
+[[fig:writer-example2-2]] since ~output_list~ is set with the updated value
+of ~fh~.
+
+#+NAME: writer-example2-img
+#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/writer-example2.png")
+  reset
+  set grid
+  set border
+  set key
+  set xtics
+  set ytics
+  set square
+  set terminal png truecolor enhanced size 640,640
+  plot '<(sed "" $PUGS_SOURCE_DIR/doc/writer-example2.gnu)' u 1:2 lw 2 t "f" w l, '<(sed "" $PUGS_SOURCE_DIR/doc/writer-example2.gnu)' u 1:3 lw 2 t "g" w l,  '<(sed "" $PUGS_SOURCE_DIR/doc/writer-example2.gnu)' u 1:4 lw 2 t "sin(pi*f)" w l
+#+END_SRC
+
+#+CAPTION: Illustration of the life time of variables. The output for ~fh~ corresponds to its value when ~output_list~ is created: the interpolation of $x -> x$.
+#+NAME: fig:writer-example2
+#+ATTR_LATEX: :width 0.38\textwidth
+#+ATTR_HTML: :width 300px;
+#+RESULTS: writer-example2-img
+
+#+NAME: writer-example2-2-img
+#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/writer-example2-2.png")
+  reset
+  set grid
+  set border
+  set key
+  set xtics
+  set ytics
+  set square
+  set terminal png truecolor enhanced size 640,640
+  plot '<(sed "" $PUGS_SOURCE_DIR/doc/writer-example2-2.gnu)' u 1:2 lw 2 t "f" w l, '<(sed "" $PUGS_SOURCE_DIR/doc/writer-example2-2.gnu)' u 1:3 lw 2 t "g" w l
+#+END_SRC
+
+#+CAPTION: Illustration of the life time of variables. ~output_list~ is updated after ~fh~ has been updated.
+#+NAME: fig:writer-example2-2
+#+ATTR_LATEX: :width 0.38\textwidth
+#+ATTR_HTML: :width 300px;
+#+RESULTS: writer-example2-2-img
+
+***** ~gnuplot~ writers <<gnuplot-writers>>
+
+There is two ~gnuplot~ writers. One is dedicated to output of dimension
+1 results (~gnuplot_1d_writer~) and the other allows post processing in
+dimension 1 and 2 (~gnuplot_writer~).
+
+Both of these writers can be used for single output or time series
+output. In the case of single output, the filename is completed by
+adding the extension ~.gnu~, in the case of time series, the filename is
+extending by adding ~.abcd.gnu~, where ~abcd~ is the number of the output
+in the series.
+
+#+BEGIN_note
+The ~gnuplot~ writers are implemented in parallel.
+
+The ~gnuplot~ post processing of produced files is the same whichever is
+the number of processors (as soon as the saved data is also the same,
+which is warrantied by ~pugs~ for explicit methods).
+#+END_note
+
+For an obvious practical reason, each ~gnuplot~ file starts with a
+preamble that indicates running information.
+
+This information contains, the production date of the file. The ~pugs~
+version (including the hash key of the last commit) and a summary of
+the output data to easy its use. For ~gnuplot~ files, this is the
+location where the provided name (to the discrete function) is used.
+
+In the case of time series the "physical" time of output is also
+stored as a comment before the variable list.
+
+Here is an example of preamble of a produced ~gnuplot~ file.
+#+NAME: head-gnuplot-txt
+#+BEGIN_SRC shell :exports results :results output
+  head -n 8 writer-example1.gnu
+#+END_SRC
+#+RESULTS: head-gnuplot-txt
+
+****** ~gnuplot_1d_writer~ functions
+
+This writer family make only sense in 1d.
+
+#+BEGIN_note
+In parallel, as soon as the saved data themselves are the same, the
+~gnuplot_1d_writer~ generates *exactly* the same files since the
+coordinates of the post processed data are sorted according to growing
+abscissa.
+#+END_note
+
+******* ~gnuplot_1d_writer: string -> writer~
+
+Defines a ~gnuplot~ writer for 1d data. This is a single file writer in
+the sense that it does not take care of time series.
+
+Let us consider an example.
+#+BEGIN_SRC pugs :exports both :results none
+  import mesh;
+  import scheme;
+  import writer;
+  import math;
+
+  let m:mesh, m = cartesianMesh([0], [1], 5);
+
+  let identity:R^1 -> R, x -> x[0];
+
+  let xh:Vh, xh = interpolate(m, P0(), identity);
+
+  let gp:writer, gp = gnuplot_1d_writer("gp_1d_example");
+  write(gp, (name_output(xh*xh, "x^2"), name_output(3*xh+2, "3x+2")));
+#+END_SRC
+
+The produced file (~"gp_1d_example.gnu"~) contains the following
+#+NAME: gp-1d-example-txt
+#+BEGIN_SRC shell :exports results :results output
+cat gp_1d_example.gnu
+#+END_SRC
+#+RESULTS: gp-1d-example-txt
+
+As one can see the output for $\mathbb{P}_0$ data consists in saving
+the values at cell centers in the order of growing abscissa. This
+allows to plot "smoother" curves than the ~gnuplot_writer~ (see
+[[gnuplot-writer-function]]).
+
+A typical use of this writer is the following.
+#+BEGIN_SRC pugs :exports both :results none
+  import mesh;
+  import scheme;
+  import writer;
+  import math;
+
+  let m:mesh, m = cartesianMesh([0], [2], 20);
+
+  let pi:R, pi = acos(-1);
+  let sin_pi_x:R^1 -> R, x -> sin(pi*x[0]);
+
+  let fh:Vh, fh = interpolate(m, P0(), sin_pi_x);
+
+  write(gnuplot_1d_writer("gp_1d_sin"), name_output(fh, "f"));
+#+END_SRC
+
+#+NAME: writer-gp-1d-sin-img
+#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/gp-1d-sin.png")
+  reset
+  set grid
+  set border
+  unset key
+  set xtics
+  set ytics
+  set square
+  set terminal png truecolor enhanced size 628,400
+  plot '<(sed "" $PUGS_SOURCE_DIR/doc/gp_1d_sin.gnu)' lw 2 w lp
+#+END_SRC
+
+#+CAPTION: Example of produced gnuplot results from the ~gnuplot_1d_writer~. Since data are stored by growing abscissa, one can use the ~w lp~ plotting option in ~gnuplot~.
+#+NAME: fig:writer-gp-1d-sin
+#+ATTR_LATEX: :width 0.38\textwidth
+#+ATTR_HTML: :width 300px;
+#+RESULTS: writer-gp-1d-sin-img
+
+
+******* ~gnuplot_1d_writer: string*R -> writer~ <<gnuplot-1d-series>>
+
+This writer differs from the previous one by handling output
+series. The real value argument defines the period to respect between
+two outputs. It can be viewed as an helper to outputs.
+
+Let us give an example to fix ideas.
+#+BEGIN_SRC pugs :exports both :results none
+  import mesh;
+  import scheme;
+  import writer;
+  import math;
+
+  let m:mesh, m = cartesianMesh([0], [1], 20);
+
+  let pi:R, pi = acos(-1);
+
+  let t:R,   t = 0;
+  let tmax:R, tmax = 1;
+  let dt:R, dt = 3.7e-2;
+
+  let period:R, period = 0.15; // sets the period output
+
+  let gp:writer, gp = gnuplot_1d_writer("gp_1d_exp_sin", period);
+
+  let f: R^1 -> R, x -> exp(-t)*sin(2*pi*x[0]);
+  let fh:Vh, fh = interpolate(m, P0(), f);
+
+  write(gp, name_output(fh, "f"), 0);
+  do {
+    if (dt > tmax-t) {
+      dt = tmax-t;
+    }
+    t +=dt;
+
+    write(gp, name_output(fh, "f"), t);
+  } while(t<tmax);
+#+END_SRC
+
+Running this example produces the following files
+#+NAME: ls-produced-gp-1d-series
+#+BEGIN_SRC shell :exports results :results output
+  ls gp_1d_exp_sin.*.gnu
+#+END_SRC
+#+RESULTS: ls-produced-gp-1d-series
+
+Each of these file contains the numerical solution at following saving
+times:
+#+NAME: times-in-gp-1d-series
+#+BEGIN_SRC shell :exports results :results output
+  grep -n "# time = " gp_1d_exp_sin.*.gnu | cut -d '=' -f 2
+#+END_SRC
+#+RESULTS: times-in-gp-1d-series
+
+****** ~gnuplot_writer~ functions <<gnuplot-writer-function>>
+
+This writer differs from the previous one since it draws the cells and
+affects the cell value to the nodes. This produces larger files but
+allow 2d representation. Also, if the saved data is exactly the same
+in parallel, the order of the cells is generally different since they
+are written processor by processor.
+
+Additionally, this writer allows to write 2d meshes, see paragraph
+[[write-mesh]].
+
+******* ~gnuplot_writer: string -> writer~
+
+Here is an illustrating example in dimension 1, see the result on
+Figure [[fig:writer-gp-sin]].
+#+BEGIN_SRC pugs :exports both :results none
+  import mesh;
+  import scheme;
+  import writer;
+  import math;
+
+  let m:mesh, m = cartesianMesh([0], [2], 20);
+
+  let pi:R, pi = acos(-1);
+  let sin_pi_x:R^1 -> R, x -> sin(pi*x[0]);
+
+  let fh:Vh, fh = interpolate(m, P0(), sin_pi_x);
+
+  write(gnuplot_writer("gp_sin"), name_output(fh, "f"));
+#+END_SRC
+
+#+NAME: writer-gp-sin-img
+#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/gp-sin.png")
+  reset
+  set grid
+  set border
+  unset key
+  set xtics
+  set ytics
+  set square
+  set terminal png truecolor enhanced size 628,400
+  plot '<(sed "" $PUGS_SOURCE_DIR/doc/gp_sin.gnu)' lw 2 w lp
+#+END_SRC
+
+#+CAPTION: Example of produced gnuplot results from the ~gnuplot_writer~. One can compare ths produced result to the one of the ~gnuplot_1d_writer~ given on Figure [[fig:writer-gp-1d-sin]]
+#+NAME: fig:writer-gp-sin
+#+ATTR_LATEX: :width 0.38\textwidth
+#+ATTR_HTML: :width 300px;
+#+RESULTS: writer-gp-sin-img
+
+Let use give a 2d example.
+#+BEGIN_SRC pugs :exports both :results none
+  import mesh;
+  import scheme;
+  import writer;
+  import math;
+
+  let m:mesh, m = cartesianMesh(-[1,1], [1,1], (40,40));
+
+  let pi:R, pi = acos(-1);
+  let f:R^2 -> R, x -> cos(pi*x[0])*sin(pi*x[1]);
+
+  write(gnuplot_writer("gp_2d_cos_sin"),
+        name_output(interpolate(m,P0(), f), "f"));
+#+END_SRC
+
+The gnuplot result is displayed on Figure [[fig:writer-gp-2d-cos-sin]].
+
+#+NAME: writer-gp-2d-cos-sin-img
+#+BEGIN_SRC gnuplot :exports results :file (substitute-in-file-name "${PUGS_SOURCE_DIR}/doc/gp-2d-cos-sin.png")
+  reset
+  set grid
+  set border
+  unset key
+  set xtics
+  set ytics
+  set square
+  unset colorbox
+  set palette rgb 33,13,10
+  set terminal png truecolor enhanced size 300,300
+  plot '<(sed "" $PUGS_SOURCE_DIR/doc/gp_2d_cos_sin.gnu)' w filledcurves closed palette, '<(sed "" $PUGS_SOURCE_DIR/doc/gp_2d_cos_sin.gnu)' lc rgb 0 w l
+#+END_SRC
+
+#+CAPTION: Example of 2d plot from the ~gnuplot_writer~
+#+NAME: fig:writer-gp-2d-cos-sin
+#+ATTR_LATEX: :width 0.38\textwidth
+#+ATTR_HTML: :width 300px;
+#+RESULTS: writer-gp-2d-cos-sin-img
+
+******* ~gnuplot_writer: string*R -> writer~
+
+This is the time series function in the case of the ~gnuplot_writer~. It
+behaves the same as [[gnuplot-1d-series]].
+
+***** ~vtk~ writers
+
+For more complex post processing (including 3d), ~pugs~ can generate ~vtk~
+outputs.
+
+The used format is one file in the ~vtu~ format for each parallel domain
+(and eventually each time). The output is done using binary data for
+performance reasons. For each time step a ~pvtu~ file is generate to
+handle parallelism. And for a complete time series, a ~pvd~ file is
+produced. This is the file that should be loaded.
+
+Observe that each of these files (~vtu~, ~pvti~ and ~pvd~) contains a
+comment that contains the creation date and the version of ~pugs~ that
+was run.
+
+The use is exactly the same than for ~gnuplot~ writers so we do not
+provide full examples.
+
+~vtk~ writers are compatible with the ~write_mesh~ function, see paragraph
+[[write-mesh]].
+
+****** ~vtk_writer: string -> writer~
+
+One should use this writer for single output (no time series).  The
+produced ~pvd~ file is built by adding ~.pvd~ to the provided ~string~.
+
+****** ~vtk_writer: string*R -> writer~
+
+This function follows the same rule. One just specifies the output
+period. The generated ~pvd~ file is built the same way, one adds ~.pvd~ to
+the provided ~string~.
+
+
+***** ~write~, ~force_write~ and ~write_mesh~ functions
+
+One a mesh writer has been defined, these functions are called to
+effectively generate the post processing files.
+
+****** ~write: writer*(output) -> void~
+
+As a parameter, it takes a ~writer~ that is not a time series one and a
+list of ~output~ (which are named discrete functions that are defined on
+the *same* mesh).
+
+There have already shown a lot of examples of use of the ~write~
+function. So let us focus on improper calls.
+
+Trying to use a time series ~writer~
+#+NAME: cannot-use-time-series-writer
+#+BEGIN_SRC pugs-error :exports both :results output
+  import mesh;
+  import scheme;
+  import writer;
+  import math;
+
+  let m:mesh, m = cartesianMesh([0], [2], 20);
+
+  let pi:R, pi = acos(-1);
+  let sin_pi_x:R^1 -> R, x -> sin(pi*x[0]);
+
+  let fh:Vh, fh = interpolate(m, P0(), sin_pi_x);
+
+  let vtk:writer, vtk = vtk_writer("vtk_time_series", 0);
+  write(vtk, name_output(fh, "f"));
+#+END_SRC
+produces the following runtime error
+#+results: cannot-use-time-series-writer
+
+Trying to use variables defined on different meshes
+#+NAME: cannot-use-diffrent-meshes-in-writer
+#+BEGIN_SRC pugs-error :exports both :results output
+  import mesh;
+  import scheme;
+  import writer;
+  import math;
+
+  let m0:mesh, m0 = cartesianMesh([0], [2], 20);
+  let m1:mesh, m1 = cartesianMesh([0], [2], 20);
+
+  let pi:R, pi = acos(-1);
+  let sin_pi_x:R^1 -> R, x -> sin(pi*x[0]);
+
+  let fh:Vh, fh = interpolate(m0, P0(), sin_pi_x);
+  let gh:Vh, gh = interpolate(m1, P0(), sin_pi_x);
+
+  let vtk:writer, vtk = vtk_writer("vtk_different_meshes");
+  write(vtk, (name_output(fh, "f"), name_output(gh, "g")));
+#+END_SRC
+gives the runtime error
+#+results: cannot-use-diffrent-meshes-in-writer
+
+****** ~write: writer*(output)*R -> void~
+
+This ~write~ function is used to save time series data. The real
+parameter (of type ~R~) is the current time.
+
+One cannot use a ~writer~ that does not support time series.
+
+Trying to use a non time series ~writer~
+#+NAME: cannot-use-non-time-series-writer
+#+BEGIN_SRC pugs-error :exports both :results output
+  import mesh;
+  import scheme;
+  import writer;
+  import math;
+
+  let m:mesh, m = cartesianMesh([0], [2], 20);
+
+  let pi:R, pi = acos(-1);
+  let sin_pi_x:R^1 -> R, x -> sin(pi*x[0]);
+
+  let fh:Vh, fh = interpolate(m, P0(), sin_pi_x);
+
+  let vtk:writer, vtk = vtk_writer("vtk_time_series");
+  write(vtk, name_output(fh, "f"), 1.2);
+#+END_SRC
+gives the runtime error
+#+results: cannot-use-non-time-series-writer
+
+****** ~force_write: writer*(output)*R -> void~
+
+One probably noticed that using the ~write~ function with a time series
+~writer~, last time of the calculation may not be written (see section
+[[gnuplot-1d-series]]). The ~force_write~ function does not check that the
+saving time has been reached. It just checks that the current time has
+not already been saved.
+
+Let us improve slightly the example given in section
+[[gnuplot-1d-series]].
+#+BEGIN_SRC pugs :exports both :results none
+  import mesh;
+  import scheme;
+  import writer;
+  import math;
+
+  let m:mesh, m = cartesianMesh([0], [1], 20);
+
+  let pi:R, pi = acos(-1);
+
+  let t:R,   t = 0;
+  let tmax:R, tmax = 1;
+  let dt:R, dt = 3.7e-2;
+
+  let period:R, period = 0.15; // sets the period output
+
+  let gp:writer, gp = gnuplot_1d_writer("gp_1d_exp_sin_force", period);
+
+  let f: R^1 -> R, x -> exp(-t)*sin(2*pi*x[0]);
+  let fh:Vh, fh = interpolate(m, P0(), f);
+
+  write(gp, name_output(fh, "f"), 0);
+  do {
+    if (dt > tmax-t) {
+      dt = tmax-t;
+    }
+    t +=dt;
+
+    fh = interpolate(m, P0(), f);
+    if (t<tmax) {
+      write(gp, name_output(fh, "f"), t);
+    } else {
+      force_write(gp, name_output(fh, "f"), t);
+    }
+  } while(t<tmax);
+#+END_SRC
+
+Running this example produces the following files
+#+NAME: ls-produced-gp-1d-series-force
+#+BEGIN_SRC shell :exports results :results output
+  ls gp_1d_exp_sin_force.*.gnu
+#+END_SRC
+#+RESULTS: ls-produced-gp-1d-series-force
+One can see the additional file.
+
+Each of these file contains the numerical solution at following saving
+times:
+#+NAME: times-in-gp-1d-series-force
+#+BEGIN_SRC shell :exports results :results output
+  grep -n "# time = " gp_1d_exp_sin_force.*.gnu | cut -d '=' -f 2
+#+END_SRC
+#+RESULTS: times-in-gp-1d-series-force
+The last post processing time is now 1.
+
+****** ~write_mesh: writer*mesh -> void~ <<write-mesh>>.
+
+This function just saves a ~mesh~ using a ~writer~. The ~gnuplot_1d_writer~
+cannot write mesh. And the ~writer~ argument must not be time series
+~writer~.
 
 [fn:pugs-def] ~pugs~: Parallel Unstructured Grid Solvers
 [fn:MPI-def] ~MPI~: Message Passing Interface
-[fn:DSL-def] ~DSL~: Domain Specific Language~
+[fn:DSL-def] ~DSL~: Domain Specific Language