diff --git a/examples/telemac3d/solit/doc/solit.tex b/examples/telemac3d/solit/doc/solit.tex
index 481c7f056ca195ffc97e9df487b673b5d2b9eda0..ad6d4236f5a366008e20fadcbf78caee5ae751aa 100644
--- a/examples/telemac3d/solit/doc/solit.tex
+++ b/examples/telemac3d/solit/doc/solit.tex
@@ -1,19 +1,13 @@
-% case name
 \chapter{Solitary wave (solit)}
-%
-% - Purpose & Description:
-%     These first two parts give reader short details about the test case,
-%     the physical phenomena involved, the geometry and specify how the numerical solution will be validated
-%
 
 \section{Purpose}
-%
+
 This test demonstrates the ability of \telemac{3D} to model the
-propagation of a solitary wave with only 2, 3 or 4 vertical levels.
+propagation of a solitary wave with only 4 vertical levels.
 In an ideal case, the wave should travel without changing its shape and amplitude.
 This study demonstrates the necessity of using the non-hydrostatic
 version of the software to simulate non-linear waves propagation.
-%
+
 \section{Notations}
 
 This test-case is the same as described in \cite{Jankowski1999}. The notations
@@ -31,27 +25,32 @@ are illustrated in Figure \ref{fig:solit_notations}.
 \section{Theory}
 
 \begin{WarningBlock}{Authorship}
-The following description is extracted from the PhD thesis of Jacek A. Jankowski \cite{Jankowski1999}.
+The following description is extracted from the PhD thesis of Jacek A. Jankowski
+\cite{Jankowski1999}.
 Only slight changes were done to the text.
 \end{WarningBlock}
 
-A solitary wave is a single elevation of the water surface above an undisturbed surrounding,
-which is neither preceded nor followed by any free surface disturbances. Neglecting
-dissipation, as well as bottom and lateral boundary shear, a solitary wave travels over a
-horizontal bottom without changing its shape and velocity.
-The accuracy of the model can be evaluated by comparing the
-amplitude and celerity of the wave with its theoretical values, as well as the deformation
+A solitary wave is a single elevation of the water surface above an undisturbed
+surrounding, which is neither preceded nor followed by any free surface
+disturbances.
+Neglecting dissipation, as well as bottom and lateral boundary shear, a solitary
+wave travels over a horizontal bottom without changing its shape and velocity.
+The accuracy of the model can be evaluated by comparing the amplitude and
+celerity of the wave with its theoretical values, as well as the deformation
 of the wave as it travels.\\
 
-There are several theories that aim at giving an approximation of the free-surface and velocity
-field for this form of non-linear finite-amplitude wave. The most known ones are the Stokes and cnoïdal
-waves theories. The latter is used here, based on Laitone \cite{Laitone1960}. According to \cite{Jankowski1999}, it
-is the most frequently used for comparative studies.
-It provides approximate formulae for the
-velocity components $u$, $w$, free surface elevation $\eta$ defined as $\eta = z + h$,
-pressure $p$ and wave celerity $c$ of a
-solitary wave with a height of $H$, on a vertical section of an infinitely long channel of an undisturbed
-depth $h$ ($z = 0$ at the free-surface, see Figure \ref{fig:solit_notations}). They read:
+There are several theories that aim at giving an approximation of the
+free-surface and velocity field for this form of non-linear finite-amplitude
+wave.
+The most known ones are the Stokes and cnoïdal waves theories.
+The latter is used here, based on Laitone \cite{Laitone1960}.
+According to \cite{Jankowski1999}, it is the most frequently used for
+comparative studies.
+It provides approximate formulae for the velocity components $u$, $w$, free
+surface elevation $\eta$ defined as $\eta = z + h$, pressure $p$ and wave
+celerity $c$ of a solitary wave with a height of $H$, on a vertical section of
+an infinitely long channel of an undisturbed depth $h$ ($z = 0$ at the
+free-surface, see Figure \ref{fig:solit_notations}). They read:
 \begin{equation}
 \label{eq:solit}
 \left\{\begin{array}{l}
@@ -63,24 +62,29 @@ p= \rho g (\eta + z) \medskip \\
 c = \sqrt{g(H+h)}
 \end{array}\right.
 \end{equation}
-It is interesting that in this analytical approximation the vertical velocity component is
-not treated as small, as commonly taken, but the pressure can be assumed hydrostatic
-(fourth line of \eqref{eq:solit}) at the same level of exactness as the horizontal velocity
-(with $\mathcal{O}((H/h)^2)$)\cite{Laitone1960}.
-Therefore, this initial condition is suitable for fair comparisons between models
-with and without hydrostatic approximation and the initial value of zero hydrodynamic
-pressure is assumed. Although the initial velocity field (first two lines of \eqref{eq:solit})
-is perfectly divergence-free, larger values of the hydrodynamic pressure appear immediately after the first time
-step (60\% of the value of the hydrostatic pressure at the bottom).\\
-
-Following the test cases provided by Ramaswamy \cite{Ramaswamy1990}, a solitary wave described by
-\eqref{eq:solit} is applied in a long channel as an initial condition, and the
-behaviour of the solution is observed thereafter. Due to the fact that the simulation is
-performed in a finite domain, and Laitone’s formulae are valid for an indefinitely long
-channel, care must be taken choosing the initial position of the wave crest. In order to
-deal with it, the use of the effective wave length $\lambda$ concept is made. $\lambda$ is equal to the
-doubled length between the wave crest and a point, where the free surface elevation is
-$\eta(x) = 0.01H$. According to Laitone:
+It is interesting that in this analytical approximation the vertical velocity
+component is not treated as small, as commonly taken, but the pressure can be
+assumed hydrostatic (fourth line of \eqref{eq:solit}) at the same level of
+exactness as the horizontal velocity (with $\mathcal{O}((H/h)^2)$)
+\cite{Laitone1960}.
+Therefore, this initial condition is suitable for fair comparisons between
+models with and without hydrostatic approximation and the initial value of zero
+hydrodynamic pressure is assumed.
+Although the initial velocity field (first two lines of \eqref{eq:solit}) is
+perfectly divergence-free, larger values of the hydrodynamic pressure appear
+immediately after the first time step (60\% of the value of the hydrostatic
+pressure at the bottom).\\
+
+Following the test cases provided by Ramaswamy \cite{Ramaswamy1990}, a solitary
+wave described by \eqref{eq:solit} is applied in a long channel as an initial
+condition, and the behaviour of the solution is observed thereafter.
+Due to the fact that the simulation is performed in a finite domain, and
+Laitone’s formulae are valid for an indefinitely long channel, care must be
+taken choosing the initial position of the wave crest.
+In order to deal with it, the use of the effective wave length $\lambda$ concept
+is made.
+$\lambda$ is equal to the doubled length between the wave crest and a point,
+where the free surface elevation is $\eta(x) = 0.01H$. According to Laitone:
 \begin{equation}
 \lambda = 6.9\sqrt{\dfrac{h^3}{H}}.
 \end{equation}
@@ -103,7 +107,8 @@ should be at least 110~m.
 \subsection{Geometry and Mesh}
 
 \subsubsection{Geometry}
-A long channel of 600~m long and 6~m wide is considered, with a constant depth $h$ = 10~m.
+A long channel of 600~m long and 6~m wide is considered, with a constant depth
+$h$ = 10~m.
 The bottom is flat, at the elevation $z$ = -10~m.
 
 \subsubsection{Mesh}
@@ -121,33 +126,42 @@ Various numbers of planes on the vertical are considered: 2, 3 and 4.
 \end{figure}
 
 \subsection{Physical parameters}
-%
-The flow is assumed to be inviscid, without shear on the lateral boundaries and on the bottom.
-%
+
+The flow is assumed to be inviscid, without shear on the lateral boundaries and
+on the bottom.
+
 \subsection{Initial and Boundary Conditions}
-%
+
 \subsubsection{Initial conditions}
-%
-As the initial condition the hydrostatic approximation
-given by \eqref{eq:solit} is applied, with a wave height of $H$ = 1~m,
-and the initial crest position at $x$ = 150~m. The velocity field and the free-surface
-elevation are given by \eqref{eq:solit}.
-%
+
+As the initial condition the hydrostatic approximation given by \eqref{eq:solit}
+is applied, with a wave height of $H$ = 1~m, and the initial crest position at
+$x$ = 150~m.
+The velocity field and the free-surface elevation are given by \eqref{eq:solit}.
+
 \subsubsection{Boundary conditions}
-%
+
 The lateral boundaries and the bottom are impervious, there is no friction.
-%
+
 \subsection{General parameters}
-%
-The time step size is constant, equal to $\Delta t$ = 0.1~s (the Courant number in the direction
-of wave propagation varies from 0.2 to about 1.0 at the wave crest).
+
+The time step size is constant, equal to $\Delta t$ = 0.1~s (the Courant number
+in the direction of wave propagation varies from 0.2 to about 1.0 at the wave
+crest).
 The simulation time is 30~s.
 
 \subsection{Numerical parameters}
-%
-The non-hydrostatic version is used. Advection of velocities is done with the method of characteristics.
-The solver for the pressure is the GMRES (Generalized Minimal Residual Method, scheme 7)
-with diagonal preconditionning, an accuracy of 10$^{-6}$ is asked for.
+
+The non-hydrostatic version is used.
+Advection of velocities is done with the method of characteristics.
+The solver for the pressure is the GMRES (Generalized Minimal Residual Method,
+scheme 7) with diagonal preconditionning, an accuracy of 10$^{-9}$ is asked for.
+
+The solver for the diffusion of velocities is the default one (conjugate
+gradient, \#1) combined with preconditioning for diffusion of velocities
+= 34 to be more efficient.
+It enables to easily reached the accuracy of 10$^{-15}$.
+As also easily reached, the accuracy for propagation is set to 10$^{-15}$.
 
 A coefficient of implicitation for the depth and velocities of 0.51 is used.
 A mass-lumping is used on the depth, with a coefficient of 1.
@@ -166,18 +180,20 @@ colours represent the velocity magnitude, from 0 (blue) to 1 (red).}
 \end{center}
 \end{figure}
 %
-Figure \ref{fig:solit_results} presents a longitudinal cross profile of the free surface
-at different times of the simulation and for the various configurations
+Figure \ref{fig:solit_results} presents a longitudinal cross profile of the free
+surface at different times of the simulation and for the various configurations
 on vertical discretisation.
 The amplitude of the wave remains nearly constant in all cases.
-The relative error on the wave amplitude after 30~s is of 4.8~\% with 2 planes and about 0.01~\%
-with 3 and 4 planes (see the Table~\ref{tab:solit_celerity}).
+The relative error on the wave amplitude after 30~s is of 4.8~\% with 2 planes,
+1~\% with 3 planes and about 0.02~\% with 4 planes
+(see the Table~\ref{tab:solit_celerity}).
 On the other hand, the theoretical celerity of the propagation is equal to $\sqrt{g(h+H)}$.
 With a water level $h$ = 10~m and a wave height $H$ = 1~m,
 the theoretical celerity is equal to 10.38~m.s$^{-1}$.
-Table~\ref{tab:solit_celerity} shows the values of the wave displacement and celerity
-on the numerical models (using 2, 3 and 4 planes). The error on the celerity is also
-displayed, showing a maximum error of 1.9~\% (using two planes) and less than 1~\% when refining on the vertical.
+Table~\ref{tab:solit_celerity} shows the values of the wave displacement and
+celerity on the numerical models (using 2, 3 and 4 planes).
+The error on the celerity is also displayed, showing a maximum error of 1.9~\%
+(using two planes) and less than 1~\% when refining on the vertical.
 \begin{table}[H]
 \caption{Solitary wave propagation over a flat bottom:
 values of the wave height and displacement after 30~s in the simulation,
@@ -188,21 +204,25 @@ using 2, 3 and 4 planes in the mesh.}
 \hline
 ~ & \textbf{2 planes} & \textbf{3 planes} & \textbf{4 planes}\\
 \hline
-\textbf{Final wave height (m)} & 0.952 & 0.999875 & 1.0002 \\
+\textbf{Final wave height (m)} & 0.9520 & 0.9900 & 1.0002 \\
+% Values read in 2D RESULT FILE for variable FREE SURFACE
 \hline
 \textbf{Final displacement (m)} & 305.5 & 309 & 310 \\
+% Find at which abscisse the crest/peak of FS is located at final time
 \hline
 \textbf{Wave celerity (m.s$^{-1}$)} & 10.18 & 10.3 & 10.33 \\
+% = Final displacement / 30 s
 \hline
 \textbf{Relative error on the celerity (\%)} & 1.9 & 0.8 & 0.5 \\
 \hline
 \end{tabular}\end{center}
 \end{table}
 \begin{WarningBlock}{Results using the hydrostatic version of TELEMAC-3D}
-The users should be aware that when using the hydrostatic version on this test-case, the
-results are not in good agreement with the approximate theoretical solution: the shape of the velocity field
-is deteriorated, the error on the wave celerity is of about 10~\% and the error on the wave height after 30~s
-of propagation is of 15~\%.
+The users should be aware that when using the hydrostatic version on this
+test-case, the results are not in good agreement with the approximate
+theoretical solution: the shape of the velocity field is deteriorated, the error
+on the wave celerity is of about 10~\% and the error on the wave height after
+30~s of propagation is of 15~\%.
 \end{WarningBlock}
 %\begin{table}[H]
 %\caption{Solitary wave propagation over a flat bottom:
@@ -227,4 +247,5 @@ of propagation is of 15~\%.
 %
 \section{Conclusion}
 %
-\telemac{3D} correctly simulates the propagation of a solitary wave when using the non-hydrostatic formulation.
+\telemac{3D} correctly simulates the propagation of a solitary wave when using
+the non-hydrostatic formulation and at least 4 planes.
diff --git a/examples/telemac3d/solit/f3d_solit_p2.slf b/examples/telemac3d/solit/f3d_solit_p2.slf
index 26fc3d98dae28b2613fbad475ab40e87b5dfa946..e1bf512b5b64367791b47a2efef426440d6e1b48 100644
--- a/examples/telemac3d/solit/f3d_solit_p2.slf
+++ b/examples/telemac3d/solit/f3d_solit_p2.slf
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:959cba383a46a679f8b145824c7ff03ebc5ce5ce81a3b62a870369eb7b95b4a3
+oid sha256:9057688b98ad1803360586bf7ccd71273643bc9201db4698e37b38e44023a399
 size 409148
diff --git a/examples/telemac3d/solit/f3d_solit_p3.slf b/examples/telemac3d/solit/f3d_solit_p3.slf
index 021e26f25a9d44b85c2c907c7aba58b2f7cba89d..72f22c82ca9fd937682db0723204c8e392ea1fe3 100644
--- a/examples/telemac3d/solit/f3d_solit_p3.slf
+++ b/examples/telemac3d/solit/f3d_solit_p3.slf
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:bb66884532f8ac5a51e7bb6bf06c93d47401b79cd738e4062a2a5cf764f6abfa
+oid sha256:0564a8baefc156f7ae754d42be1767cbd5419011c841add8e9354f0aec961ce2
 size 699972
diff --git a/examples/telemac3d/solit/f3d_solit_p4.slf b/examples/telemac3d/solit/f3d_solit_p4.slf
index 2cb46b70ed2b71385348c0b9a052ed9d71bb8519..e4f179cd855cf2b42808365c5ea1948494c0108b 100644
--- a/examples/telemac3d/solit/f3d_solit_p4.slf
+++ b/examples/telemac3d/solit/f3d_solit_p4.slf
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:d4fa01fba9549c749029887999d95b18753188f5033a586abfa0fd9d1b7e12a3
+oid sha256:7627510a1ad76bb4b9b3b86c9218b51f0e8a9f6effb9f94aead68eb834ccec6a
 size 990796
diff --git a/examples/telemac3d/solit/t3d_solit_p2.cas b/examples/telemac3d/solit/t3d_solit_p2.cas
index 097ca0d3f5cf0d00074d558a2d9853f8c9199258..afc8427141ff115263d24cf549ed23b6ae833c23 100644
--- a/examples/telemac3d/solit/t3d_solit_p2.cas
+++ b/examples/telemac3d/solit/t3d_solit_p2.cas
@@ -39,9 +39,9 @@ INITIAL CONDITIONS = 'SPECIAL'
 /
 SCHEME FOR ADVECTION OF VELOCITIES = 1
 /
-SOLVER FOR DIFFUSION OF VELOCITIES          = 7
-/
-ACCURACY FOR DIFFUSION OF VELOCITIES = 1.E-6
+/SOLVER FOR DIFFUSION OF VELOCITIES          = 7
+ACCURACY FOR DIFFUSION OF VELOCITIES = 1.E-15
+PRECONDITIONING FOR DIFFUSION OF VELOCITIES = 34
 /
 IMPLICITATION FOR DEPTH = 0.51
 IMPLICITATION FOR VELOCITIES = 0.51
@@ -54,8 +54,10 @@ COEFFICIENT FOR VERTICAL DIFFUSION OF VELOCITIES   : 0.
 MASS-BALANCE : YES
 INFORMATION ABOUT MASS-BALANCE FOR EACH LISTING PRINTOUT : NO
 /
+ACCURACY FOR PROPAGATION = 1.E-15
+/
 /PRECONDITIONING FOR PPE = 7  SOLVER FOR PPE = 6
-ACCURACY FOR PPE = 1.E-6
+ACCURACY FOR PPE = 1.E-9
 /
 MASS-LUMPING FOR DEPTH : 1.
 /
diff --git a/examples/telemac3d/solit/t3d_solit_p3.cas b/examples/telemac3d/solit/t3d_solit_p3.cas
index 22dcb99e776ef3a7b678d1b8abcb190760f956c5..5d60ea26533ac4c01103b1d84d05ac0566850e35 100644
--- a/examples/telemac3d/solit/t3d_solit_p3.cas
+++ b/examples/telemac3d/solit/t3d_solit_p3.cas
@@ -32,9 +32,9 @@ INITIAL CONDITIONS = 'SPECIAL'
 /
 SCHEME FOR ADVECTION OF VELOCITIES = 1
 /
-SOLVER FOR DIFFUSION OF VELOCITIES          = 7
-/
-ACCURACY FOR DIFFUSION OF VELOCITIES = 1.E-6
+/SOLVER FOR DIFFUSION OF VELOCITIES          = 7
+ACCURACY FOR DIFFUSION OF VELOCITIES = 1.E-15
+PRECONDITIONING FOR DIFFUSION OF VELOCITIES = 34
 /
 IMPLICITATION FOR DEPTH = 0.51
 IMPLICITATION FOR VELOCITIES = 0.51
@@ -47,8 +47,10 @@ COEFFICIENT FOR VERTICAL DIFFUSION OF VELOCITIES   : 0.
 MASS-BALANCE : YES
 INFORMATION ABOUT MASS-BALANCE FOR EACH LISTING PRINTOUT : NO
 /
+ACCURACY FOR PROPAGATION = 1.E-15
+/
 /PRECONDITIONING FOR PPE = 7  SOLVER FOR PPE = 6
-ACCURACY FOR PPE = 1.E-6
+ACCURACY FOR PPE = 1.E-9
 /
 MASS-LUMPING FOR DEPTH : 1.
 /
diff --git a/examples/telemac3d/solit/t3d_solit_p4.cas b/examples/telemac3d/solit/t3d_solit_p4.cas
index 578567dded728e1d38e49f245efa320bc62af30a..9d9f716bf6629b5d94cf6da53caa2284c29067ac 100644
--- a/examples/telemac3d/solit/t3d_solit_p4.cas
+++ b/examples/telemac3d/solit/t3d_solit_p4.cas
@@ -31,9 +31,9 @@ INITIAL CONDITIONS = 'SPECIAL'
 /
 SCHEME FOR ADVECTION OF VELOCITIES = 1
 /
-SOLVER FOR DIFFUSION OF VELOCITIES          = 7
-/
-ACCURACY FOR DIFFUSION OF VELOCITIES = 1.E-6
+/SOLVER FOR DIFFUSION OF VELOCITIES          = 7
+ACCURACY FOR DIFFUSION OF VELOCITIES = 1.E-15
+PRECONDITIONING FOR DIFFUSION OF VELOCITIES = 34
 /
 IMPLICITATION FOR DEPTH = 0.51
 IMPLICITATION FOR VELOCITIES = 0.51
@@ -46,8 +46,10 @@ COEFFICIENT FOR VERTICAL DIFFUSION OF VELOCITIES   : 0.
 MASS-BALANCE : YES
 INFORMATION ABOUT MASS-BALANCE FOR EACH LISTING PRINTOUT : NO
 /
+ACCURACY FOR PROPAGATION = 1.E-15
+/
 /PRECONDITIONING FOR PPE = 7  SOLVER FOR PPE = 6
-ACCURACY FOR PPE = 1.E-6
+ACCURACY FOR PPE = 1.E-9
 /
 MASS-LUMPING FOR DEPTH : 1.
 /
diff --git a/examples/telemac3d/solit/vnv_solit.py b/examples/telemac3d/solit/vnv_solit.py
index 5dd563cfd742f360196378fa6a4dbb843bad7011..2e6f1161f625a1e006f8ddbb92477bb60ee2c980 100644
--- a/examples/telemac3d/solit/vnv_solit.py
+++ b/examples/telemac3d/solit/vnv_solit.py
@@ -86,32 +86,32 @@ class VnvStudy(AbstractVnvStudy):
         # Comparison with the last time frame of the reference file.
         self.check_epsilons('vnv_1:T3DRES',
                             'f3d_solit_p2.slf',
-                            eps=[1.E-10])
+                            eps=[1.E-12])
 
         # Comparison with the last time frame of the reference file.
         self.check_epsilons('vnv_2:T3DRES',
                             'f3d_solit_p2.slf',
-                            eps=[1.E-10])
+                            eps=[1.E-12])
 
         # Comparison between sequential and parallel run.
         self.check_epsilons('vnv_1:T3DRES',
                             'vnv_2:T3DRES',
-                            eps=[1.E-10])
+                            eps=[1.E-12])
 
         # Comparison with the last time frame of the reference file.
         self.check_epsilons('vnv_3:T3DRES',
                             'f3d_solit_p3.slf',
-                            eps=[1.E-10])
+                            eps=[1.E-11])
 
         # Comparison with the last time frame of the reference file.
         self.check_epsilons('vnv_4:T3DRES',
                             'f3d_solit_p3.slf',
-                            eps=[1.E-10])
+                            eps=[1.E-11])
 
         # Comparison between sequential and parallel run.
         self.check_epsilons('vnv_3:T3DRES',
                             'vnv_4:T3DRES',
-                            eps=[1.E-10])
+                            eps=[1.E-11])
 
         # Comparison with the last time frame of the reference file.
         self.check_epsilons('vnv_5:T3DRES',
@@ -151,6 +151,7 @@ class VnvStudy(AbstractVnvStudy):
                    filled_contours=True,
                    xlim=(1, 600),
                    ylim=(-10, 0),
+                   y_label='z (m)',
                    fig_size=(20, 4),
                    fig_name='img/solit_2planes')
 
@@ -163,6 +164,7 @@ class VnvStudy(AbstractVnvStudy):
                    filled_contours=True,
                    xlim=(1, 600),
                    ylim=(-10, 0),
+                   y_label='z (m)',
                    fig_size=(20, 4),
                    fig_name='img/solit_3planes')
 
@@ -175,6 +177,7 @@ class VnvStudy(AbstractVnvStudy):
                    filled_contours=True,
                    xlim=(1, 600),
                    ylim=(-10, 0),
+                   y_label='z (m)',
                    fig_size=(20, 4),
                    fig_name='img/solit_4planes')