diff --git a/test/douar_stdout.txt b/test/douar_stdout.txt
new file mode 100644
index 0000000000000000000000000000000000000000..37f4fa19a4d337b8726da7026d3143c218de03a0
--- /dev/null
+++ b/test/douar_stdout.txt
@@ -0,0 +1,1523 @@
+ ----------------------------------------------------------------------- 
+ ----------------------------------------------------------------------- 
+ program called with no argument
+ ----------------------------------------------------------------------- 
+ ----------------------------------------------------------------------- 
+               Start of Computations  Total time   Step time   Increment
+                       Reading Input      0.0002      0.0002      0.0002
+                       Problem Setup      0.2935      0.2935      0.2934
+                     define surfaces      0.2937      0.2937      0.0001
+                          surface 01      0.2937      0.2937      0.0000
+                        define cloud      0.3095      0.3095      0.0158
+              define velocity octree      0.3097      0.3097      0.0002
+ ----------------------------------------------------------------------- 
+ ----------------------------------------------------------------------- 
+               Start of Step       1      0.3195
+ ----------------------------------------------------------------------- 
+ ----------------------------------------------------------------------- 
+                  Loop over surfaces      0.3197      0.0002      0.0002
+                          surface 01      0.3197      0.0002      0.0000
+                                                                        S. 1:     4356points
+                      refine surface      0.3197      0.0003      0.0000
+                                                                        S. 1:   0 added ptcls in refine_surface
+                      check delaunay      0.3337      0.0142      0.0139
+                                                                        S. 1:    0 flips in delaunay2
+                                                                        S. 1:     4356points
+ ----------------------------------------------------------------------- 
+          Start of iteration       1      0.3379
+ ----------------------------------------------------------------------- 
+                 Create octree solve      0.3380      0.0002      0.0002
+                                                                        (1) params%griditer < 0                F
+                                                                        (2) current_level==params%levelmax_oct F
+                                                                        (3) increase_current_level             T
+                                                                        converge = (1) & (2) & (3) ->  F
+        compute ref. criterion on ov      0.3396      0.0018      0.0016
+                                                                        osolve%nleaves=    4096
+        improve osolve for criterion      0.3411      0.0032      0.0015
+                                                                        max(crit)        =   0.000000E+00
+                                                                        max(crit,w_l_i_f)=      -Infinity
+                                                                        osolve%nleaves=    4096
+    Improve osolve for imposed boxes      0.3412      0.0034      0.0001
+                                                                        osolve%nleaves=  434176
+    Smoothen octree after refinement      0.6582      0.3203      0.3170
+                                                                        osolve%nleaves=  434176
+                   Build osolve icon      2.1942      1.8563      1.5360
+             Imbed surface in osolve      4.0485      3.7107      1.8543
+                embedding surface  1      4.0488      3.7109      0.0002
+ ----------------------------------------------------------------------- 
+ ----------------------------------------------------------------------- 
+ program called with no argument
+ ----------------------------------------------------------------------- 
+ ----------------------------------------------------------------------- 
+               Start of Computations  Total time   Step time   Increment
+                       Reading Input      0.0001      0.0001      0.0001
+                       Problem Setup      0.3187      0.3187      0.3186
+                     define surfaces      0.3187      0.3187      0.0000
+                          surface 01      0.3187      0.3187      0.0000
+                        define cloud      0.3318      0.3318      0.0131
+              define velocity octree      0.3319      0.3319      0.0001
+ ----------------------------------------------------------------------- 
+ ----------------------------------------------------------------------- 
+               Start of Step       1      0.3416
+ ----------------------------------------------------------------------- 
+ ----------------------------------------------------------------------- 
+                  Loop over surfaces      0.3417      0.0001      0.0001
+                          surface 01      0.3417      0.0001      0.0000
+                                                                        S. 1:     4356points
+                      refine surface      0.3417      0.0001      0.0000
+                                                                        S. 1:   0 added ptcls in refine_surface
+                      check delaunay      0.3468      0.0051      0.0051
+                                                                        S. 1:    0 flips in delaunay2
+                                                                        S. 1:     4356points
+ ----------------------------------------------------------------------- 
+          Start of iteration       1      0.3512
+ ----------------------------------------------------------------------- 
+                 Create octree solve      0.3513      0.0001      0.0001
+                                                                        (1) params%griditer < 0                F
+                                                                        (2) current_level==params%levelmax_oct F
+                                                                        (3) increase_current_level             T
+                                                                        converge = (1) & (2) & (3) ->  F
+        compute ref. criterion on ov      0.3528      0.0016      0.0016
+                                                                        osolve%nleaves=    4096
+        improve osolve for criterion      0.3544      0.0032      0.0016
+                                                                        max(crit)        =   0.000000E+00
+                                                                        max(crit,w_l_i_f)=      -Infinity
+                                                                        osolve%nleaves=    4096
+    Improve osolve for imposed boxes      0.3545      0.0033      0.0001
+                                                                        osolve%nleaves=  434176
+    Smoothen octree after refinement      0.6716      0.3204      0.3171
+                                                                        osolve%nleaves=  434176
+                   Build osolve icon      2.1915      1.8403      1.5200
+             Imbed surface in osolve      3.8104      3.4591      1.6188
+                embedding surface  1      3.8104      3.4592      0.0001
+                                                                        osolve%nleaves=  434176
+          Assessing octree stability    132.8824    132.5312    129.0720
+                                                                        current refine level     4
+                                                                        After criterion based refinement:
+                                                                        nleaves_old=        0 leaves
+                                                                        nleaves_new=   434176 leaves
+                                                                        After surfaces embedding:
+                                                                        nleaves_old=        0 leaves
+                                                                        nleaves_new=   434176 leaves
+                                                                        C2: authorise increase of refine level =  F
+                     Change 3D cloud    132.9096    132.5584      0.0272
+        Interpolate velo onto osolve    149.2604    148.9092     16.3507
+                                                                        ov%pressure:      0.0000E+00 0.0000E+00
+                                                                        osolve%pressure:  0.0000E+00 0.0000E+00
+                      Find bad faces    160.5220    160.1708     11.2617
+                           Find void    162.4849    162.1337      1.9629
+                                                                        whole leaf in fluid   196608
+          Define boundary conditions    162.5888    162.2376      0.1039
+                         wsmp setup1    162.6166    162.2654      0.0279
+                                                                        n =   698922
+                                                                        nz= 27029961
+                         wsmp setup2    167.7975    167.4463      5.1809
+ -----------------------------------
+ start of non-linear iteratio      1    168.0293
+ -----------------------------------
+                        build system    168.0294      0.0001      0.0001
+                                                                        nelem per proc (min/max)       27059       42683
+                                                                        viscosity range   1.00000E-04  3.59904E+02
+                                                                        viscosity capped in    212992
+                          wsmp solve    215.1212     47.0919     47.0918
+                                                                        =================================
+                                                                        wsmp: initialisation time:         0.3170919 s
+                                                                        wsmp: ordering time:              20.2699749 s
+                                                                        wsmp: symbolic fact. time:         4.0266092 s
+                                                                        wsmp: Cholesky fact. time:       165.7119761 s
+                                                                        wsmp: Factorisation            16063.1933253 Mflops
+                                                                        wsmp: back substitution time:      0.9624178 s
+                                                                        wsmp: from b to rhs vector:        0.0383561 s
+                                                                        ---------------------------------
+                                                                        wsmp: total solve time:          191.3325701 s
+                                                                        =================================
+                                                                        u : -0.80899E-01  0.11421E+01
+                                                                        v : -0.12881E+00  0.18421E+00
+                                                                        w : -0.38890E+00  0.45830E+00
+                                                                        =================================
+                 Calculate pressures    406.5988    238.5694    191.4775
+                                                                        raw    pressures : -0.15257E+04  0.15258E+04
+                 Smoothing pressures    406.9503    238.9210      0.3515
+                                                                        smooth pressures : -0.37450E+03  0.30944E+01
+                do leaf measurements    407.0367    239.0074      0.0864
+       compute convergence criterion    407.3943    239.3650      0.3576
+                                                                        velocity_diff_norm = 0.5710589 (tol = 0.0100000)
+                                                                        increase_current_level -> F
+                                                                        current_level -> 4
+                    slicing the cube    407.4115    239.3822      0.0172
+Compute isostasy and adjust vertical    407.4122    239.3829      0.0006
+                  Compute divergence    407.4122    239.3829      0.0000
+                        wsmp_cleanup    407.6406    239.6113      0.2284
+              Reset surface geometry    407.6507    239.6213      0.0101
+ -----------------------------------------------------------------------
+           Temperature calculations     407.6514    239.6221      0.0007
+ -----------------------------------------------------------------------
+                         wsmp setup1    407.6718    239.6425      0.0204
+                                                                        n =   232974
+                                                                        nz=  5929000
+                         wsmp setup2    408.7394    240.7101      1.0676
+                        build system    408.8176    240.7883      0.0782
+                                                                        nelem per proc (min/max)       27059       42683
+                                                                        ===================================
+                                                                        wsmp: initialisation time:         0.0000811 s
+                                                                        wsmp: analysis time:               7.5014892 s
+                                                                        wsmp: nb of exp. non-z in factors   182869452
+                                                                        wsmp: nb of exp. Mflops in factorisation 207504.4011855
+                                                                        wsmp: LU time:                    16.8921990 s
+                                                                        wsmp: actual nb of non-z in factors   182869452
+                                                                        wsmp: actual nb of Mflops in factorisation  12287.4867499
+                                                                        wsmp: back substitution time:      0.2088301 s
+                                                                        wsmp: from b to rhs vector:        0.0121360 s
+                                                                        ---------------------------------
+                                                                        wsmp: total solve time:           24.8069861 s
+                                                                        =================================
+                                                                        temp :  0.00000E+00  0.11187E+01
+                                                                        =================================
+                        wsmp_cleanup    440.4120    272.3827     31.5944
+                 Advect the surfaces    440.6456    272.6162      0.2335
+                   Advect surface  1    440.6456    272.6163      0.0001
+                                                                        S. 1:    7 flips in delaunay2
+                                                                        S. 1:    0 flips in delaunay2
+                    Erode surface  1    440.7042    272.6749      0.0586
+                 Update cloud fields    469.9109    301.8816     29.2068
+                    Advect the cloud    475.2267    307.1973      5.3157
+                        Write output    506.7627    338.7334     31.5360
+                    End of time step    517.9713    349.9420     11.2086
+ ----------------------------------------------------------------------- 
+ ----------------------------------------------------------------------- 
+               Start of Step       2    517.9714
+ ----------------------------------------------------------------------- 
+ ----------------------------------------------------------------------- 
+                  Loop over surfaces    517.9714      0.0000      0.0000
+                          surface 01    517.9714      0.0000      0.0000
+                                                                        S. 1:     4356points
+                      refine surface    517.9714      0.0001      0.0000
+                                                                        S. 1:   0 added ptcls in refine_surface
+                      check delaunay    517.9811      0.0097      0.0096
+                                                                        S. 1:   35 flips in delaunay2
+                                                                        S. 1:    0 flips in delaunay2
+                                                                        S. 1:     4356points
+ ----------------------------------------------------------------------- 
+          Start of iteration       1    517.9860
+ ----------------------------------------------------------------------- 
+                 Create octree solve    517.9861      0.0001      0.0001
+                                                                        (1) params%griditer < 0                F
+                                                                        (2) current_level==params%levelmax_oct F
+                                                                        (3) increase_current_level             F
+                                                                        converge = (1) & (2) & (3) ->  F
+        compute ref. criterion on ov    517.9877      0.0016      0.0016
+                                                                        osolve%nleaves=    4096
+        improve osolve for criterion    518.1387      0.1527      0.1511
+                                                                        max(crit)        =   6.660066E+01
+                                                                        max(crit,w_l_i_f)=   6.527429E+01
+                                                                        osolve%nleaves=    4096
+    Improve osolve for imposed boxes    518.2569      0.2709      0.1182
+                                                                        osolve%nleaves=  434176
+    Smoothen octree after refinement    518.5738      0.5878      0.3170
+                                                                        osolve%nleaves=  434176
+                   Build osolve icon    520.0938      2.1078      1.5200
+             Imbed surface in osolve    521.8237      3.8377      1.7299
+                embedding surface  1    521.8238      3.8378      0.0001
+                                                                        osolve%nleaves=  434176
+          Assessing octree stability    650.8906    132.9046    129.0669
+                                                                        current refine level     4
+                                                                        After criterion based refinement:
+                                                                        nleaves_old=   434176 leaves
+                                                                        nleaves_new=   434176 leaves
+                                                                        After surfaces embedding:
+                                                                        nleaves_old=   434176 leaves
+                                                                        nleaves_new=   434176 leaves
+                                                                        C2: authorise increase of refine level =  T
+                     Change 3D cloud    650.9142    132.9282      0.0235
+        Interpolate velo onto osolve    667.8279    149.8419     16.9137
+                                                                        ov%pressure:     -0.5053E+03 0.4650E+01
+                                                                        osolve%pressure: -0.3745E+03 0.3094E+01
+                      Find bad faces    678.7312    160.7452     10.9034
+                           Find void    680.6563    162.6703      1.9251
+                                                                        whole leaf in fluid   196608
+          Define boundary conditions    680.7564    162.7704      0.1001
+                         wsmp setup1    680.7863    162.8003      0.0300
+                                                                        n =   698922
+                                                                        nz= 27029961
+                         wsmp setup2    685.9745    167.9885      5.1882
+ -----------------------------------
+ start of non-linear iteratio      1    686.2039
+ -----------------------------------
+                        build system    686.2039      0.0001      0.0001
+                                                                        nelem per proc (min/max)       27059       42683
+                                                                        viscosity range   1.00000E-04  4.73035E+01
+                                                                        viscosity capped in    212992
+                          wsmp solve    734.4003     48.1964     48.1963
+                                                                        =================================
+                                                                        wsmp: initialisation time:         0.2316258 s
+                                                                        wsmp: ordering time:              20.3056488 s
+                                                                        wsmp: symbolic fact. time:         4.0109699 s
+                                                                        wsmp: Cholesky fact. time:       165.3999331 s
+                                                                        wsmp: Factorisation            16093.4980913 Mflops
+                                                                        wsmp: back substitution time:      0.9578860 s
+                                                                        wsmp: from b to rhs vector:        0.0371420 s
+                                                                        ---------------------------------
+                                                                        wsmp: total solve time:          190.9443481 s
+                                                                        =================================
+                                                                        u : -0.90323E-02  0.11161E+01
+                                                                        v : -0.14402E+00  0.17181E+00
+                                                                        w : -0.47633E+00  0.50369E+00
+                                                                        =================================
+                 Calculate pressures    925.4953    239.2915    191.0951
+                                                                        raw    pressures : -0.11579E+03  0.83714E+02
+                 Smoothing pressures    925.8496    239.6457      0.3543
+                                                                        smooth pressures : -0.65802E+02  0.24927E+00
+                do leaf measurements    925.9531    239.7492      0.1034
+       compute convergence criterion    926.3165    240.1126      0.3634
+                                                                        velocity_diff_norm = 0.1672622 (tol = 0.0100000)
+                                                                        increase_current_level -> T
+                                                                        current_level -> 4
+                    slicing the cube    926.3337    240.1298      0.0172
+Compute isostasy and adjust vertical    926.3344    240.1306      0.0007
+                  Compute divergence    926.3345    240.1306      0.0000
+                        wsmp_cleanup    926.5624    240.3585      0.2279
+              Reset surface geometry    926.5745    240.3706      0.0121
+ -----------------------------------------------------------------------
+           Temperature calculations     926.5751    240.3712      0.0006
+ -----------------------------------------------------------------------
+                         wsmp setup1    926.5958    240.3919      0.0207
+                                                                        n =   232974
+                                                                        nz=  5929000
+                         wsmp setup2    927.6460    241.4421      1.0501
+                        build system    927.7205    241.5166      0.0746
+                                                                        nelem per proc (min/max)       27059       42683
+                                                                        ===================================
+                                                                        wsmp: initialisation time:         0.0000830 s
+                                                                        wsmp: analysis time:               7.4842472 s
+                                                                        wsmp: nb of exp. non-z in factors   182869452
+                                                                        wsmp: nb of exp. Mflops in factorisation 207504.4011855
+                                                                        wsmp: LU time:                    16.9641008 s
+                                                                        wsmp: actual nb of non-z in factors   182869452
+                                                                        wsmp: actual nb of Mflops in factorisation  12235.4066307
+                                                                        wsmp: back substitution time:      0.2020431 s
+                                                                        wsmp: from b to rhs vector:        0.0121160 s
+                                                                        ---------------------------------
+                                                                        wsmp: total solve time:           24.9322319 s
+                                                                        =================================
+                                                                        temp :  0.00000E+00  0.12351E+01
+                                                                        =================================
+                        wsmp_cleanup    959.5596    273.3557     31.8390
+                 Advect the surfaces    959.7915    273.5876      0.2319
+                   Advect surface  1    959.7915    273.5876      0.0000
+                                                                        S. 1:    4 flips in delaunay2
+                                                                        S. 1:    0 flips in delaunay2
+                    Erode surface  1    959.8403    273.6364      0.0488
+                 Update cloud fields    983.8068    297.6029     23.9665
+                    Advect the cloud    989.3985    303.1946      5.5917
+                        Write output   1022.1014    335.8975     32.7029
+                    End of time step   1033.5868    347.3829     11.4854
+ ----------------------------------------------------------------------- 
+ ----------------------------------------------------------------------- 
+               Start of Step       3   1033.5869
+ ----------------------------------------------------------------------- 
+ ----------------------------------------------------------------------- 
+                  Loop over surfaces   1033.5869      0.0000      0.0000
+                          surface 01   1033.5869      0.0000      0.0000
+                                                                        S. 1:     4356points
+                      refine surface   1033.5869      0.0001      0.0000
+                                                                        S. 1:   0 added ptcls in refine_surface
+                      check delaunay   1033.5944      0.0076      0.0075
+                                                                        S. 1:    5 flips in delaunay2
+                                                                        S. 1:    0 flips in delaunay2
+                                                                        S. 1:     4356points
+ ----------------------------------------------------------------------- 
+          Start of iteration       1   1033.5989
+ ----------------------------------------------------------------------- 
+                 Create octree solve   1033.5990      0.0001      0.0001
+                                                                        (1) params%griditer < 0                F
+                                                                        (2) current_level==params%levelmax_oct F
+                                                                        (3) increase_current_level             T
+                                                                        converge = (1) & (2) & (3) ->  F
+        compute ref. criterion on ov   1033.6006      0.0016      0.0016
+                                                                        osolve%nleaves=    4096
+        improve osolve for criterion   1033.7482      0.1493      0.1476
+                                                                        max(crit)        =   6.628113E+01
+                                                                        max(crit,w_l_i_f)=   5.610013E+01
+                                                                        osolve%nleaves=    4096
+    Improve osolve for imposed boxes   1033.8664      0.2675      0.1183
+                                                                        osolve%nleaves=  434176
+    Smoothen octree after refinement   1034.1824      0.5835      0.3160
+                                                                        osolve%nleaves=  434176
+                   Build osolve icon   1035.7056      2.1067      1.5232
+             Imbed surface in osolve   1037.4169      3.8180      1.7112
+                embedding surface  1   1037.4170      3.8180      0.0001
+                                                                        osolve%nleaves=  434176
+          Assessing octree stability   1166.8541    133.2552    129.4372
+                                                                        current refine level     4
+                                                                        After criterion based refinement:
+                                                                        nleaves_old=   434176 leaves
+                                                                        nleaves_new=   434176 leaves
+                                                                        After surfaces embedding:
+                                                                        nleaves_old=   434176 leaves
+                                                                        nleaves_new=   434176 leaves
+                                                                        C2: authorise increase of refine level =  T
+                     Change 3D cloud   1166.8788    133.2799      0.0247
+        Interpolate velo onto osolve   1184.2274    150.6284     17.3485
+                                                                        ov%pressure:     -0.8490E+02 0.4586E+00
+                                                                        osolve%pressure: -0.6580E+02 0.2493E+00
+                      Find bad faces   1191.6915    158.0926      7.4642
+                           Find void   1193.6233    160.0244      1.9318
+                                                                        whole leaf in fluid   196608
+          Define boundary conditions   1193.7252    160.1263      0.1019
+                         wsmp setup1   1193.7575    160.1586      0.0323
+                                                                        n =   698922
+                                                                        nz= 27029961
+                         wsmp setup2   1198.9579    165.3590      5.2004
+ -----------------------------------
+ start of non-linear iteratio      1   1199.1861
+ -----------------------------------
+                        build system   1199.1861      0.0001      0.0001
+                                                                        nelem per proc (min/max)       27059       42683
+                                                                        viscosity range   1.00000E-04  8.65793E+00
+                                                                        viscosity capped in    212992
+                          wsmp solve   1246.7021     47.5160     47.5160
+                                                                        =================================
+                                                                        wsmp: initialisation time:         0.3375721 s
+                                                                        wsmp: ordering time:              20.2010119 s
+                                                                        wsmp: symbolic fact. time:         4.0342140 s
+                                                                        wsmp: Cholesky fact. time:       164.2253020 s
+                                                                        wsmp: Factorisation            16208.6077817 Mflops
+                                                                        wsmp: back substitution time:      0.9580081 s
+                                                                        wsmp: from b to rhs vector:        0.0368600 s
+                                                                        ---------------------------------
+                                                                        wsmp: total solve time:          189.7942071 s
+                                                                        =================================
+                                                                        u :  0.00000E+00  0.10831E+01
+                                                                        v : -0.12618E+00  0.26206E+00
+                                                                        w : -0.40763E+00  0.52970E+00
+                                                                        =================================
+                 Calculate pressures   1436.6457    237.4596    189.9436
+                                                                        raw    pressures : -0.54550E+02  0.54404E+02
+                 Smoothing pressures   1437.0126    237.8265      0.3669
+                                                                        smooth pressures : -0.14630E+02  0.00000E+00
+                do leaf measurements   1437.1094    237.9233      0.0968
+       compute convergence criterion   1437.4683    238.2823      0.3589
+                                                                        velocity_diff_norm = 0.1075870 (tol = 0.0100000)
+                                                                        increase_current_level -> T
+                                                                        current_level -> 4
+                    slicing the cube   1437.4857    238.2996      0.0174
+Compute isostasy and adjust vertical   1437.4866    238.3005      0.0009
+                  Compute divergence   1437.4866    238.3005      0.0000
+                        wsmp_cleanup   1437.7183    238.5322      0.2317
+              Reset surface geometry   1437.7309    238.5449      0.0126
+ -----------------------------------------------------------------------
+           Temperature calculations    1437.7316    238.5455      0.0006
+ -----------------------------------------------------------------------
+                         wsmp setup1   1437.7520    238.5659      0.0204
+                                                                        n =   232974
+                                                                        nz=  5929000
+                         wsmp setup2   1438.8030    239.6169      1.0510
+                        build system   1438.8782    239.6921      0.0751
+                                                                        nelem per proc (min/max)       27059       42683
+                                                                        ===================================
+                                                                        wsmp: initialisation time:         0.0000811 s
+                                                                        wsmp: analysis time:               7.4918990 s
+                                                                        wsmp: nb of exp. non-z in factors   182869452
+                                                                        wsmp: nb of exp. Mflops in factorisation 207504.4011855
+                                                                        wsmp: LU time:                    16.9326000 s
+                                                                        wsmp: actual nb of non-z in factors   182869452
+                                                                        wsmp: actual nb of Mflops in factorisation  12258.1689529
+                                                                        wsmp: back substitution time:      0.2025290 s
+                                                                        wsmp: from b to rhs vector:        0.0120542 s
+                                                                        ---------------------------------
+                                                                        wsmp: total solve time:           24.9643869 s
+                                                                        =================================
+                                                                        temp :  0.00000E+00  0.13135E+01
+                                                                        =================================
+                        wsmp_cleanup   1470.6924    271.5063     31.8143
+                 Advect the surfaces   1470.9260    271.7400      0.2336
+                   Advect surface  1   1470.9261    271.7400      0.0000
+                                                                        S. 1:    1 flips in delaunay2
+                                                                        S. 1:    0 flips in delaunay2
+                    Erode surface  1   1470.9747    271.7887      0.0486
+                 Update cloud fields   1493.3929    294.2068     22.4182
+                    Advect the cloud   1499.0517    299.8656      5.6588
+                        Write output   1532.3497    333.1636     33.2980
+                    End of time step   1544.1431    344.9570     11.7934
+ ----------------------------------------------------------------------- 
+ ----------------------------------------------------------------------- 
+               Start of Step       4   1544.1431
+ ----------------------------------------------------------------------- 
+ ----------------------------------------------------------------------- 
+                  Loop over surfaces   1544.1432      0.0000      0.0000
+                          surface 01   1544.1432      0.0000      0.0000
+                                                                        S. 1:     4356points
+                      refine surface   1544.1432      0.0001      0.0000
+                                                                        S. 1:   0 added ptcls in refine_surface
+                      check delaunay   1544.1509      0.0078      0.0077
+                                                                        S. 1:    1 flips in delaunay2
+                                                                        S. 1:    0 flips in delaunay2
+                                                                        S. 1:     4356points
+ ----------------------------------------------------------------------- 
+          Start of iteration       1   1544.1554
+ ----------------------------------------------------------------------- 
+                 Create octree solve   1544.1555      0.0001      0.0001
+                                                                        (1) params%griditer < 0                F
+                                                                        (2) current_level==params%levelmax_oct F
+                                                                        (3) increase_current_level             T
+                                                                        converge = (1) & (2) & (3) ->  F
+        compute ref. criterion on ov   1544.1571      0.0017      0.0016
+                                                                        osolve%nleaves=    4096
+        improve osolve for criterion   1544.3065      0.1511      0.1494
+                                                                        max(crit)        =   6.512954E+01
+                                                                        max(crit,w_l_i_f)=   5.792272E+01
+                                                                        osolve%nleaves=    4096
+    Improve osolve for imposed boxes   1544.4329      0.2775      0.1264
+                                                                        osolve%nleaves=  434176
+    Smoothen octree after refinement   1544.7497      0.5942      0.3168
+                                                                        osolve%nleaves=  434176
+                   Build osolve icon   1546.2738      2.1184      1.5241
+             Imbed surface in osolve   1547.9820      3.8265      1.7082
+                embedding surface  1   1547.9820      3.8266      0.0001
+                                                                        osolve%nleaves=  434176
+          Assessing octree stability   1677.4144    133.2590    129.4324
+                                                                        current refine level     4
+                                                                        After criterion based refinement:
+                                                                        nleaves_old=   434176 leaves
+                                                                        nleaves_new=   434176 leaves
+                                                                        After surfaces embedding:
+                                                                        nleaves_old=   434176 leaves
+                                                                        nleaves_new=   434176 leaves
+                                                                        C2: authorise increase of refine level =  T
+                     Change 3D cloud   1677.4345    133.2791      0.0201
+        Interpolate velo onto osolve   1694.9915    150.8361     17.5570
+                                                                        ov%pressure:     -0.1676E+02 0.7618E-01
+                                                                        osolve%pressure: -0.1463E+02 0.0000E+00
+                      Find bad faces   1700.0431    155.8877      5.0516
+                           Find void   1701.9784    157.8230      1.9353
+                                                                        whole leaf in fluid   196608
+          Define boundary conditions   1702.0763    157.9209      0.0979
+                         wsmp setup1   1702.1023    157.9469      0.0260
+                                                                        n =   698922
+                                                                        nz= 27029961
+                         wsmp setup2   1707.2777    163.1223      5.1754
+ -----------------------------------
+ start of non-linear iteratio      1   1707.5054
+ -----------------------------------
+                        build system   1707.5054      0.0000      0.0000
+                                                                        nelem per proc (min/max)       27059       42683
+                                                                        viscosity range   2.25826E-04  2.77222E+00
+                                                                        viscosity capped in    212992
+                          wsmp solve   1755.0529     47.5475     47.5475
+                                                                        =================================
+                                                                        wsmp: initialisation time:         0.3607311 s
+                                                                        wsmp: ordering time:              20.3277299 s
+                                                                        wsmp: symbolic fact. time:         4.0358078 s
+                                                                        wsmp: Cholesky fact. time:       164.3769121 s
+                                                                        wsmp: Factorisation            16193.6580591 Mflops
+                                                                        wsmp: back substitution time:      0.9633529 s
+                                                                        wsmp: from b to rhs vector:        0.0370040 s
+                                                                        ---------------------------------
+                                                                        wsmp: total solve time:          190.1028478 s
+                                                                        =================================
+                                                                        u : -0.47155E-02  0.10807E+01
+                                                                        v : -0.12969E+00  0.19265E+00
+                                                                        w : -0.35192E+00  0.53567E+00
+                                                                        =================================
+                 Calculate pressures   1945.3068    237.8015    190.2540
+                                                                        raw    pressures : -0.42853E+02  0.42715E+02
+                 Smoothing pressures   1945.6632    238.1578      0.3564
+                                                                        smooth pressures : -0.32806E+01  0.00000E+00
+                do leaf measurements   1945.7586    238.2532      0.0954
+       compute convergence criterion   1946.1078    238.6024      0.3492
+                                                                        velocity_diff_norm = 0.0730174 (tol = 0.0100000)
+                                                                        increase_current_level -> T
+                                                                        current_level -> 4
+                    slicing the cube   1946.1251    238.6198      0.0173
+Compute isostasy and adjust vertical   1946.1261    238.6207      0.0009
+                  Compute divergence   1946.1261    238.6208      0.0001
+                        wsmp_cleanup   1946.3639    238.8585      0.2377
+              Reset surface geometry   1946.3754    238.8701      0.0116
+ -----------------------------------------------------------------------
+           Temperature calculations    1946.3763    238.8709      0.0008
+ -----------------------------------------------------------------------
+                         wsmp setup1   1946.3995    238.8942      0.0233
+                                                                        n =   232974
+                                                                        nz=  5929000
+                         wsmp setup2   1947.4511    239.9458      1.0516
+                        build system   1947.5301    240.0247      0.0789
+                                                                        nelem per proc (min/max)       27059       42683
+                                                                        ===================================
+                                                                        wsmp: initialisation time:         0.0000820 s
+                                                                        wsmp: analysis time:               7.5116599 s
+                                                                        wsmp: nb of exp. non-z in factors   182869452
+                                                                        wsmp: nb of exp. Mflops in factorisation 207504.4011855
+                                                                        wsmp: LU time:                    16.9493830 s
+                                                                        wsmp: actual nb of non-z in factors   182869452
+                                                                        wsmp: actual nb of Mflops in factorisation  12246.0311166
+                                                                        wsmp: back substitution time:      0.2026129 s
+                                                                        wsmp: from b to rhs vector:        0.0124860 s
+                                                                        ---------------------------------
+                                                                        wsmp: total solve time:           25.0780189 s
+                                                                        =================================
+                                                                        temp :  0.00000E+00  0.14017E+01
+                                                                        =================================
+                        wsmp_cleanup   1979.5815    272.0761     32.0514
+                 Advect the surfaces   1979.8160    272.3106      0.2345
+                   Advect surface  1   1979.8162    272.3108      0.0001
+                                                                        S. 1:    0 flips in delaunay2
+                    Erode surface  1   1979.8637    272.3583      0.0475
+                 Update cloud fields   1998.8362    291.3308     18.9726
+                    Advect the cloud   2004.6513    297.1459      5.8151
+                        Write output   2038.4955    330.9901     33.8442
+                    End of time step   2050.2551    342.7497     11.7596
+ ----------------------------------------------------------------------- 
+ ----------------------------------------------------------------------- 
+               Start of Step       5   2050.2553
+ ----------------------------------------------------------------------- 
+ ----------------------------------------------------------------------- 
+                  Loop over surfaces   2050.2554      0.0001      0.0001
+                          surface 01   2050.2554      0.0001      0.0000
+                                                                        S. 1:     4356points
+                      refine surface   2050.2555      0.0001      0.0000
+                                                                        S. 1:   0 added ptcls in refine_surface
+                      check delaunay   2050.2620      0.0067      0.0066
+                                                                        S. 1:    0 flips in delaunay2
+                                                                        S. 1:     4356points
+ ----------------------------------------------------------------------- 
+          Start of iteration       1   2050.2653
+ ----------------------------------------------------------------------- 
+                 Create octree solve   2050.2654      0.0002      0.0002
+                                                                        (1) params%griditer < 0                F
+                                                                        (2) current_level==params%levelmax_oct F
+                                                                        (3) increase_current_level             T
+                                                                        converge = (1) & (2) & (3) ->  F
+        compute ref. criterion on ov   2050.2670      0.0018      0.0016
+                                                                        osolve%nleaves=    4096
+        improve osolve for criterion   2050.4125      0.1472      0.1455
+                                                                        max(crit)        =   6.436618E+01
+                                                                        max(crit,w_l_i_f)=   5.140083E+01
+                                                                        osolve%nleaves=    4096
+    Improve osolve for imposed boxes   2050.5374      0.2722      0.1249
+                                                                        osolve%nleaves=  434176
+    Smoothen octree after refinement   2050.8503      0.5850      0.3129
+                                                                        osolve%nleaves=  434176
+                   Build osolve icon   2052.3739      2.1086      1.5236
+             Imbed surface in osolve   2054.1735      3.9083      1.7997
+                embedding surface  1   2054.1737      3.9085      0.0002
+                                                                        osolve%nleaves=  434176
+          Assessing octree stability   2183.8168    133.5515    129.6431
+                                                                        current refine level     4
+                                                                        After criterion based refinement:
+                                                                        nleaves_old=   434176 leaves
+                                                                        nleaves_new=   434176 leaves
+                                                                        After surfaces embedding:
+                                                                        nleaves_old=   434176 leaves
+                                                                        nleaves_new=   434176 leaves
+                                                                        C2: authorise increase of refine level =  T
+                     Change 3D cloud   2183.8304    133.5651      0.0136
+        Interpolate velo onto osolve   2201.5362    151.2709     17.7058
+                                                                        ov%pressure:     -0.3487E+01 0.0000E+00
+                                                                        osolve%pressure: -0.3281E+01 0.0000E+00
+                      Find bad faces   2202.7813    152.5161      1.2451
+                           Find void   2204.7295    154.4643      1.9482
+                                                                        whole leaf in fluid   196608
+          Define boundary conditions   2204.8316    154.5663      0.1020
+                         wsmp setup1   2204.8611    154.5959      0.0296
+                                                                        n =   698922
+                                                                        nz= 27029961
+                         wsmp setup2   2210.0533    159.7880      5.1921
+ -----------------------------------
+ start of non-linear iteratio      1   2210.2913
+ -----------------------------------
+                        build system   2210.2914      0.0001      0.0001
+                                                                        nelem per proc (min/max)       27059       42683
+                                                                        viscosity range   3.09806E-04  1.26972E+00
+                                                                        viscosity capped in    212992
+                          wsmp solve   2257.5288     47.2375     47.2374
+                                                                        =================================
+                                                                        wsmp: initialisation time:         0.4169340 s
+                                                                        wsmp: ordering time:              20.1170051 s
+                                                                        wsmp: symbolic fact. time:         3.9956121 s
+                                                                        wsmp: Cholesky fact. time:       163.2416360 s
+                                                                        wsmp: Factorisation            16306.2780565 Mflops
+                                                                        wsmp: back substitution time:      0.9628792 s
+                                                                        wsmp: from b to rhs vector:        0.0369651 s
+                                                                        ---------------------------------
+                                                                        wsmp: total solve time:          188.7722960 s
+                                                                        =================================
+                                                                        u : -0.11237E-01  0.10543E+01
+                                                                        v : -0.13328E+00  0.14841E+00
+                                                                        w : -0.27844E+00  0.50564E+00
+                                                                        =================================
+                 Calculate pressures   2446.4491    236.1578    188.9203
+                                                                        raw    pressures : -0.25934E+02  0.25784E+02
+                 Smoothing pressures   2446.8133    236.5220      0.3642
+                                                                        smooth pressures : -0.79245E+00  0.00000E+00
+                do leaf measurements   2446.9029    236.6116      0.0895
+       compute convergence criterion   2447.2734    236.9821      0.3705
+                                                                        velocity_diff_norm = 0.0594301 (tol = 0.0100000)
+                                                                        increase_current_level -> T
+                                                                        current_level -> 4
+                    slicing the cube   2447.2906    236.9993      0.0172
+Compute isostasy and adjust vertical   2447.2914    237.0001      0.0008
+                  Compute divergence   2447.2915    237.0002      0.0001
+                        wsmp_cleanup   2447.5199    237.2287      0.2284
+              Reset surface geometry   2447.5315    237.2402      0.0116
+ -----------------------------------------------------------------------
+           Temperature calculations    2447.5324    237.2411      0.0008
+ -----------------------------------------------------------------------
+                         wsmp setup1   2447.5522    237.2609      0.0199
+                                                                        n =   232974
+                                                                        nz=  5929000
+                         wsmp setup2   2448.6019    238.3106      1.0497
+                        build system   2448.6702    238.3789      0.0682
+                                                                        nelem per proc (min/max)       27059       42683
+                                                                        ===================================
+                                                                        wsmp: initialisation time:         0.0000801 s
+                                                                        wsmp: analysis time:               7.4665198 s
+                                                                        wsmp: nb of exp. non-z in factors   182869452
+                                                                        wsmp: nb of exp. Mflops in factorisation 207504.4011855
+                                                                        wsmp: LU time:                    16.9795361 s
+                                                                        wsmp: actual nb of non-z in factors   182869452
+                                                                        wsmp: actual nb of Mflops in factorisation  12224.2840548
+                                                                        wsmp: back substitution time:      0.1998639 s
+                                                                        wsmp: from b to rhs vector:        0.0122809 s
+                                                                        ---------------------------------
+                                                                        wsmp: total solve time:           25.1174109 s
+                                                                        =================================
+                                                                        temp :  0.00000E+00  0.14788E+01
+                                                                        =================================
+                        wsmp_cleanup   2480.6938    270.4025     32.0236
+                 Advect the surfaces   2480.9274    270.6361      0.2336
+                   Advect surface  1   2480.9275    270.6363      0.0002
+                                                                        S. 1:    1 flips in delaunay2
+                                                                        S. 1:    0 flips in delaunay2
+                    Erode surface  1   2480.9760    270.6847      0.0485
+                 Update cloud fields   2498.4937    288.2024     17.5177
+                    Advect the cloud   2504.3671    294.0758      5.8734
+                        Write output   2538.5407    328.2494     34.1735
+                    End of time step   2550.4081    340.1168     11.8674
+ ----------------------------------------------------------------------- 
+ ----------------------------------------------------------------------- 
+               Start of Step       6   2550.4083
+ ----------------------------------------------------------------------- 
+ ----------------------------------------------------------------------- 
+                  Loop over surfaces   2550.4084      0.0001      0.0001
+                          surface 01   2550.4084      0.0001      0.0000
+                                                                        S. 1:     4356points
+                      refine surface   2550.4084      0.0001      0.0000
+                                                                        S. 1:   0 added ptcls in refine_surface
+                      check delaunay   2550.4139      0.0056      0.0054
+                                                                        S. 1:    0 flips in delaunay2
+                                                                        S. 1:     4356points
+ ----------------------------------------------------------------------- 
+          Start of iteration       1   2550.4172
+ ----------------------------------------------------------------------- 
+                 Create octree solve   2550.4174      0.0002      0.0002
+                                                                        (1) params%griditer < 0                F
+                                                                        (2) current_level==params%levelmax_oct F
+                                                                        (3) increase_current_level             T
+                                                                        converge = (1) & (2) & (3) ->  F
+        compute ref. criterion on ov   2550.4190      0.0018      0.0016
+                                                                        osolve%nleaves=    4096
+        improve osolve for criterion   2550.5647      0.1474      0.1456
+                                                                        max(crit)        =   6.409702E+01
+                                                                        max(crit,w_l_i_f)=   5.137951E+01
+                                                                        osolve%nleaves=    4096
+    Improve osolve for imposed boxes   2550.6935      0.2763      0.1289
+                                                                        osolve%nleaves=  434176
+    Smoothen octree after refinement   2551.0096      0.5924      0.3161
+                                                                        osolve%nleaves=  434176
+                   Build osolve icon   2552.5483      2.1311      1.5387
+             Imbed surface in osolve   2554.3220      3.9048      1.7738
+                embedding surface  1   2554.3222      3.9050      0.0002
+                                                                        osolve%nleaves=  434176
+          Assessing octree stability   2734.1773    183.7600    179.8550
+                                                                        current refine level     4
+                                                                        After criterion based refinement:
+                                                                        nleaves_old=   434176 leaves
+                                                                        nleaves_new=   434176 leaves
+                                                                        After surfaces embedding:
+                                                                        nleaves_old=   434176 leaves
+                                                                        nleaves_new=   434176 leaves
+                                                                        C2: authorise increase of refine level =  T
+                     Change 3D cloud   2734.1932    183.7760      0.0160
+        Interpolate velo onto osolve   2749.4678    199.0505     15.2745
+                                                                        ov%pressure:     -0.8091E+00 0.0000E+00
+                                                                        osolve%pressure: -0.7924E+00 0.0000E+00
+                      Find bad faces   2749.9240    199.5067      0.4562
+                           Find void   2751.8836    201.4664      1.9596
+                                                                        whole leaf in fluid   196608
+          Define boundary conditions   2751.9908    201.5736      0.1072
+                         wsmp setup1   2752.0191    201.6019      0.0283
+                                                                        n =   698922
+                                                                        nz= 27029961
+                         wsmp setup2   2757.1903    206.7730      5.1712
+ -----------------------------------
+ start of non-linear iteratio      1   2757.4212
+ -----------------------------------
+                        build system   2757.4214      0.0002      0.0002
+                                                                        nelem per proc (min/max)       27059       42683
+                                                                        viscosity range   2.98698E-04  9.49640E-01
+                                                                        viscosity capped in    212992
+                          wsmp solve   2805.6284     48.2072     48.2070
+                                                                        =================================
+                                                                        wsmp: initialisation time:         0.5335970 s
+                                                                        wsmp: ordering time:              20.1718328 s
+                                                                        wsmp: symbolic fact. time:         4.0157630 s
+                                                                        wsmp: Cholesky fact. time:       163.8344381 s
+                                                                        wsmp: Factorisation            16247.2770606 Mflops
+                                                                        wsmp: back substitution time:      0.9706621 s
+                                                                        wsmp: from b to rhs vector:        0.0368619 s
+                                                                        ---------------------------------
+                                                                        wsmp: total solve time:          189.5643780 s
+                                                                        =================================
+                                                                        u : -0.76373E-02  0.10269E+01
+                                                                        v : -0.12324E+00  0.13864E+00
+                                                                        w : -0.21613E+00  0.44693E+00
+                                                                        =================================
+                 Calculate pressures   2995.3417    237.9205    189.7133
+                                                                        raw    pressures : -0.12629E+02  0.12495E+02
+                 Smoothing pressures   2995.7043    238.2831      0.3625
+                                                                        smooth pressures : -0.23767E+00  0.00000E+00
+                do leaf measurements   2995.8056    238.3844      0.1014
+       compute convergence criterion   2996.1546    238.7334      0.3489
+                                                                        velocity_diff_norm = 0.0646290 (tol = 0.0100000)
+                                                                        increase_current_level -> T
+                                                                        current_level -> 4
+                    slicing the cube   2996.1719    238.7507      0.0173
+Compute isostasy and adjust vertical   2996.1729    238.7517      0.0010
+                  Compute divergence   2996.1729    238.7517      0.0001
+                        wsmp_cleanup   2996.4027    238.9815      0.2298
+              Reset surface geometry   2996.4151    238.9939      0.0124
+ -----------------------------------------------------------------------
+           Temperature calculations    2996.4159    238.9947      0.0008
+ -----------------------------------------------------------------------
+                         wsmp setup1   2996.4364    239.0152      0.0205
+                                                                        n =   232974
+                                                                        nz=  5929000
+                         wsmp setup2   2997.4836    240.0624      1.0472
+                        build system   2997.5509    240.1297      0.0673
+                                                                        nelem per proc (min/max)       27059       42683
+                                                                        ===================================
+                                                                        wsmp: initialisation time:         0.0000958 s
+                                                                        wsmp: analysis time:               7.4147830 s
+                                                                        wsmp: nb of exp. non-z in factors   182869452
+                                                                        wsmp: nb of exp. Mflops in factorisation 207504.4011855
+                                                                        wsmp: LU time:                    17.0323141 s
+                                                                        wsmp: actual nb of non-z in factors   182869452
+                                                                        wsmp: actual nb of Mflops in factorisation  12186.4046845
+                                                                        wsmp: back substitution time:      0.2059109 s
+                                                                        wsmp: from b to rhs vector:        0.0123830 s
+                                                                        ---------------------------------
+                                                                        wsmp: total solve time:           25.1437089 s
+                                                                        =================================
+                                                                        temp :  0.00000E+00  0.15469E+01
+                                                                        =================================
+                        wsmp_cleanup   3029.5116    272.0904     31.9607
+                 Advect the surfaces   3029.7442    272.3230      0.2325
+                   Advect surface  1   3029.7443    272.3232      0.0002
+                                                                        S. 1:    4 flips in delaunay2
+                                                                        S. 1:    0 flips in delaunay2
+                    Erode surface  1   3029.7930    272.3718      0.0486
+                 Update cloud fields   3047.7421    290.3209     17.9491
+                    Advect the cloud   3053.6112    296.1900      5.8691
+                        Write output   3088.3393    330.9181     34.7281
+                    End of time step   3100.2905    342.8693     11.9511
+ ----------------------------------------------------------------------- 
+ ----------------------------------------------------------------------- 
+               Start of Step       7   3100.2907
+ ----------------------------------------------------------------------- 
+ ----------------------------------------------------------------------- 
+                  Loop over surfaces   3100.2908      0.0001      0.0001
+                          surface 01   3100.2908      0.0001      0.0000
+                                                                        S. 1:     4356points
+                      refine surface   3100.2909      0.0002      0.0001
+                                                                        S. 1:   0 added ptcls in refine_surface
+                      check delaunay   3100.2980      0.0073      0.0071
+                                                                        S. 1:    1 flips in delaunay2
+                                                                        S. 1:    0 flips in delaunay2
+                                                                        S. 1:     4356points
+ ----------------------------------------------------------------------- 
+          Start of iteration       1   3100.3025
+ ----------------------------------------------------------------------- 
+                 Create octree solve   3100.3026      0.0002      0.0002
+                                                                        (1) params%griditer < 0                F
+                                                                        (2) current_level==params%levelmax_oct F
+                                                                        (3) increase_current_level             T
+                                                                        converge = (1) & (2) & (3) ->  F
+        compute ref. criterion on ov   3100.3042      0.0018      0.0016
+                                                                        osolve%nleaves=    4096
+        improve osolve for criterion   3100.4437      0.1412      0.1394
+                                                                        max(crit)        =   6.413314E+01
+                                                                        max(crit,w_l_i_f)=   5.057416E+01
+                                                                        osolve%nleaves=    4096
+    Improve osolve for imposed boxes   3100.5675      0.2651      0.1239
+                                                                        osolve%nleaves=  434176
+    Smoothen octree after refinement   3100.8804      0.5779      0.3128
+                                                                        osolve%nleaves=  434176
+                   Build osolve icon   3102.4034      2.1010      1.5231
+             Imbed surface in osolve   3104.1684      3.8660      1.7650
+                embedding surface  1   3104.1687      3.8662      0.0002
+                                                                        osolve%nleaves=  434176
+          Assessing octree stability   3313.2467    212.9443    209.0781
+                                                                        current refine level     4
+                                                                        After criterion based refinement:
+                                                                        nleaves_old=   434176 leaves
+                                                                        nleaves_new=   434176 leaves
+                                                                        After surfaces embedding:
+                                                                        nleaves_old=   434176 leaves
+                                                                        nleaves_new=   434176 leaves
+                                                                        C2: authorise increase of refine level =  T
+                     Change 3D cloud   3313.2627    212.9602      0.0159
+        Interpolate velo onto osolve   3331.3349    231.0324     18.0722
+                                                                        ov%pressure:     -0.2545E+00 0.0000E+00
+                                                                        osolve%pressure: -0.2377E+00 0.0000E+00
+                      Find bad faces   3331.7853    231.4829      0.4504
+                           Find void   3333.7413    233.4388      1.9560
+                                                                        whole leaf in fluid   196608
+          Define boundary conditions   3333.8375    233.5350      0.0962
+                         wsmp setup1   3333.8668    233.5643      0.0293
+                                                                        n =   698922
+                                                                        nz= 27029961
+                         wsmp setup2   3339.0627    238.7602      5.1959
+ -----------------------------------
+ start of non-linear iteratio      1   3339.2946
+ -----------------------------------
+                        build system   3339.2948      0.0002      0.0002
+                                                                        nelem per proc (min/max)       27059       42683
+                                                                        viscosity range   3.16015E-04  1.39427E+00
+                                                                        viscosity capped in    212992
+                          wsmp solve   3387.3099     48.0153     48.0151
+                                                                        =================================
+                                                                        wsmp: initialisation time:         0.5982368 s
+                                                                        wsmp: ordering time:              20.1400878 s
+                                                                        wsmp: symbolic fact. time:         4.0331759 s
+                                                                        wsmp: Cholesky fact. time:       164.0228140 s
+                                                                        wsmp: Factorisation            16228.6174840 Mflops
+                                                                        wsmp: back substitution time:      1.2452710 s
+                                                                        wsmp: from b to rhs vector:        0.0368302 s
+                                                                        ---------------------------------
+                                                                        wsmp: total solve time:          190.0774920 s
+                                                                        =================================
+                                                                        u : -0.57582E-03  0.10049E+01
+                                                                        v : -0.12076E+00  0.12935E+00
+                                                                        w : -0.18961E+00  0.40227E+00
+                                                                        =================================
+                 Calculate pressures   3577.5770    238.2823    190.2671
+                                                                        raw    pressures : -0.67642E+01  0.66260E+01
+                 Smoothing pressures   3577.9223    238.6277      0.3453
+                                                                        smooth pressures : -0.19466E+00  0.00000E+00
+                do leaf measurements   3578.0176    238.7230      0.0953
+       compute convergence criterion   3578.3694    239.0748      0.3518
+                                                                        velocity_diff_norm = 0.0742502 (tol = 0.0100000)
+                                                                        increase_current_level -> T
+                                                                        current_level -> 4
+                    slicing the cube   3578.3864    239.0918      0.0170
+Compute isostasy and adjust vertical   3578.3873    239.0926      0.0009
+                  Compute divergence   3578.3873    239.0927      0.0001
+                        wsmp_cleanup   3578.6164    239.3218      0.2291
+              Reset surface geometry   3578.6274    239.3327      0.0110
+ -----------------------------------------------------------------------
+           Temperature calculations    3578.6282    239.3336      0.0009
+ -----------------------------------------------------------------------
+                         wsmp setup1   3578.6487    239.3541      0.0205
+                                                                        n =   232974
+                                                                        nz=  5929000
+                         wsmp setup2   3579.6975    240.4029      1.0488
+                        build system   3579.7701    240.4755      0.0726
+                                                                        nelem per proc (min/max)       27059       42683
+                                                                        ===================================
+                                                                        wsmp: initialisation time:         0.0000851 s
+                                                                        wsmp: analysis time:               7.5326600 s
+                                                                        wsmp: nb of exp. non-z in factors   182869452
+                                                                        wsmp: nb of exp. Mflops in factorisation 207504.4011855
+                                                                        wsmp: LU time:                    17.0063190 s
+                                                                        wsmp: actual nb of non-z in factors   182869452
+                                                                        wsmp: actual nb of Mflops in factorisation  12205.0322185
+                                                                        wsmp: back substitution time:      0.2020421 s
+                                                                        wsmp: from b to rhs vector:        0.0123019 s
+                                                                        ---------------------------------
+                                                                        wsmp: total solve time:           25.3288782 s
+                                                                        =================================
+                                                                        temp :  0.00000E+00  0.16081E+01
+                                                                        =================================
+                        wsmp_cleanup   3611.9315    272.6369     32.1614
+                 Advect the surfaces   3612.1638    272.8692      0.2322
+                   Advect surface  1   3612.1639    272.8693      0.0002
+                                                                        S. 1:    4 flips in delaunay2
+                                                                        S. 1:    0 flips in delaunay2
+                    Erode surface  1   3612.2125    272.9178      0.0485
+                 Update cloud fields   3626.7526    287.4580     14.5402
+                    Advect the cloud   3632.6651    293.3705      5.9125
+                        Write output   3667.5310    328.2364     34.8659
+                    End of time step   3679.6790    340.3844     12.1480
+ ----------------------------------------------------------------------- 
+ ----------------------------------------------------------------------- 
+               Start of Step       8   3679.6792
+ ----------------------------------------------------------------------- 
+ ----------------------------------------------------------------------- 
+                  Loop over surfaces   3679.6793      0.0001      0.0001
+                          surface 01   3679.6793      0.0001      0.0000
+                                                                        S. 1:     4356points
+                      refine surface   3679.6794      0.0001      0.0000
+                                                                        S. 1:   0 added ptcls in refine_surface
+                      check delaunay   3679.6861      0.0068      0.0067
+                                                                        S. 1:    2 flips in delaunay2
+                                                                        S. 1:    0 flips in delaunay2
+                                                                        S. 1:     4356points
+ ----------------------------------------------------------------------- 
+          Start of iteration       1   3679.6907
+ ----------------------------------------------------------------------- 
+                 Create octree solve   3679.6908      0.0002      0.0002
+                                                                        (1) params%griditer < 0                F
+                                                                        (2) current_level==params%levelmax_oct F
+                                                                        (3) increase_current_level             T
+                                                                        converge = (1) & (2) & (3) ->  F
+        compute ref. criterion on ov   3679.6925      0.0018      0.0016
+                                                                        osolve%nleaves=    4096
+        improve osolve for criterion   3679.8320      0.1413      0.1395
+                                                                        max(crit)        =   6.406980E+01
+                                                                        max(crit,w_l_i_f)=   4.854481E+01
+                                                                        osolve%nleaves=    4096
+    Improve osolve for imposed boxes   3679.9567      0.2660      0.1247
+                                                                        osolve%nleaves=  434176
+    Smoothen octree after refinement   3680.2729      0.5822      0.3162
+                                                                        osolve%nleaves=  434176
+                   Build osolve icon   3681.8094      2.1187      1.5365
+             Imbed surface in osolve   3683.5693      3.8786      1.7599
+                embedding surface  1   3683.5695      3.8788      0.0002
+                                                                        osolve%nleaves=  434176
+          Assessing octree stability   3879.8927    200.2021    196.3232
+                                                                        current refine level     4
+                                                                        After criterion based refinement:
+                                                                        nleaves_old=   434176 leaves
+                                                                        nleaves_new=   434176 leaves
+                                                                        After surfaces embedding:
+                                                                        nleaves_old=   434176 leaves
+                                                                        nleaves_new=   434176 leaves
+                                                                        C2: authorise increase of refine level =  T
+                     Change 3D cloud   3879.9114    200.2207      0.0186
+        Interpolate velo onto osolve   3895.6023    215.9116     15.6909
+                                                                        ov%pressure:     -0.2287E+00 0.0000E+00
+                                                                        osolve%pressure: -0.1947E+00 0.0000E+00
+                      Find bad faces   3901.8023    222.1117      6.2001
+                           Find void   3903.7893    224.0986      1.9869
+                                                                        whole leaf in fluid   196608
+          Define boundary conditions   3903.8920    224.2013      0.1027
+                         wsmp setup1   3903.9308    224.2402      0.0389
+                                                                        n =   698922
+                                                                        nz= 27029961
+                         wsmp setup2   3909.1026    229.4120      5.1718
+ -----------------------------------
+ start of non-linear iteratio      1   3909.3353
+ -----------------------------------
+                        build system   3909.3354      0.0002      0.0002
+                                                                        nelem per proc (min/max)       27059       42683
+                                                                        viscosity range   3.27029E-04  7.11221E+00
+                                                                        viscosity capped in    212992
+                          wsmp solve   3956.8223     47.4870     47.4868
+                                                                        =================================
+                                                                        wsmp: initialisation time:         0.6620622 s
+                                                                        wsmp: ordering time:              20.9286959 s
+                                                                        wsmp: symbolic fact. time:         4.0489881 s
+                                                                        wsmp: Cholesky fact. time:       183.4286239 s
+                                                                        wsmp: Factorisation            14511.7127896 Mflops
+                                                                        wsmp: back substitution time:      4.3394611 s
+                                                                        wsmp: from b to rhs vector:        0.0370772 s
+                                                                        ---------------------------------
+                                                                        wsmp: total solve time:          213.4477439 s
+                                                                        =================================
+                                                                        u :  0.00000E+00  0.10004E+01
+                                                                        v : -0.11743E+00  0.11994E+00
+                                                                        w : -0.17484E+00  0.39201E+00
+                                                                        =================================
+                 Calculate pressures   4170.5535    261.2182    213.7312
+                                                                        raw    pressures : -0.45989E+01  0.44571E+01
+                 Smoothing pressures   4170.8877    261.5524      0.3342
+                                                                        smooth pressures : -0.17973E+00  0.00000E+00
+                do leaf measurements   4170.9749    261.6397      0.0872
+       compute convergence criterion   4171.3273    261.9920      0.3523
+                                                                        velocity_diff_norm = 0.0555724 (tol = 0.0100000)
+                                                                        increase_current_level -> T
+                                                                        current_level -> 4
+                    slicing the cube   4171.3447    262.0095      0.0174
+Compute isostasy and adjust vertical   4171.3455    262.0102      0.0008
+                  Compute divergence   4171.3455    262.0103      0.0001
+                        wsmp_cleanup   4171.5768    262.2415      0.2312
+              Reset surface geometry   4171.5884    262.2531      0.0116
+ -----------------------------------------------------------------------
+           Temperature calculations    4171.5892    262.2540      0.0008
+ -----------------------------------------------------------------------
+                         wsmp setup1   4171.6091    262.2739      0.0199
+                                                                        n =   232974
+                                                                        nz=  5929000
+                         wsmp setup2   4172.6590    263.3238      1.0499
+                        build system   4172.7350    263.3997      0.0759
+                                                                        nelem per proc (min/max)       27059       42683
+                                                                        ===================================
+                                                                        wsmp: initialisation time:         0.0000861 s
+                                                                        wsmp: analysis time:               7.5584779 s
+                                                                        wsmp: nb of exp. non-z in factors   182869452
+                                                                        wsmp: nb of exp. Mflops in factorisation 207504.4011855
+                                                                        wsmp: LU time:                    17.0183580 s
+                                                                        wsmp: actual nb of non-z in factors   182869452
+                                                                        wsmp: actual nb of Mflops in factorisation  12196.3982641
+                                                                        wsmp: back substitution time:      0.2110701 s
+                                                                        wsmp: from b to rhs vector:        0.0123980 s
+                                                                        ---------------------------------
+                                                                        wsmp: total solve time:           25.5622339 s
+                                                                        =================================
+                                                                        temp :  0.00000E+00  0.16637E+01
+                                                                        =================================
+                        wsmp_cleanup   4205.1416    295.8063     32.4066
+                 Advect the surfaces   4205.3771    296.0419      0.2356
+                   Advect surface  1   4205.3773    296.0420      0.0002
+                                                                        S. 1:    1 flips in delaunay2
+                                                                        S. 1:    0 flips in delaunay2
+                    Erode surface  1   4205.4258    296.0906      0.0485
+                 Update cloud fields   4219.9153    310.5800     14.4894
+                    Advect the cloud   4225.9529    316.6177      6.0377
+                        Write output   4261.0746    351.7393     35.1216
+                    End of time step   4273.1999    363.8647     12.1254
+ ----------------------------------------------------------------------- 
+ ----------------------------------------------------------------------- 
+               Start of Step       9   4273.2001
+ ----------------------------------------------------------------------- 
+ ----------------------------------------------------------------------- 
+                  Loop over surfaces   4273.2002      0.0001      0.0001
+                          surface 01   4273.2002      0.0001      0.0000
+                                                                        S. 1:     4356points
+                      refine surface   4273.2003      0.0002      0.0001
+                                                                        S. 1:   0 added ptcls in refine_surface
+                      check delaunay   4273.2076      0.0075      0.0073
+                                                                        S. 1:    1 flips in delaunay2
+                                                                        S. 1:    0 flips in delaunay2
+                                                                        S. 1:     4356points
+ ----------------------------------------------------------------------- 
+          Start of iteration       1   4273.2119
+ ----------------------------------------------------------------------- 
+                 Create octree solve   4273.2121      0.0002      0.0002
+                                                                        (1) params%griditer < 0                F
+                                                                        (2) current_level==params%levelmax_oct F
+                                                                        (3) increase_current_level             T
+                                                                        converge = (1) & (2) & (3) ->  F
+        compute ref. criterion on ov   4273.2137      0.0018      0.0016
+                                                                        osolve%nleaves=    4096
+        improve osolve for criterion   4273.3556      0.1437      0.1419
+                                                                        max(crit)        =   6.398082E+01
+                                                                        max(crit,w_l_i_f)=   4.664698E+01
+                                                                        osolve%nleaves=    4096
+    Improve osolve for imposed boxes   4273.4775      0.2656      0.1219
+                                                                        osolve%nleaves=  434176
+    Smoothen octree after refinement   4273.7906      0.5787      0.3130
+                                                                        osolve%nleaves=  434176
+                   Build osolve icon   4275.3152      2.1033      1.5246
+             Imbed surface in osolve   4277.0969      3.8850      1.7817
+                embedding surface  1   4277.0972      3.8853      0.0003
+                                                                        osolve%nleaves=  434176
+          Assessing octree stability   4447.2282    174.0163    170.1310
+                                                                        current refine level     4
+                                                                        After criterion based refinement:
+                                                                        nleaves_old=   434176 leaves
+                                                                        nleaves_new=   434176 leaves
+                                                                        After surfaces embedding:
+                                                                        nleaves_old=   434176 leaves
+                                                                        nleaves_new=   434176 leaves
+                                                                        C2: authorise increase of refine level =  T
+                     Change 3D cloud   4447.2463    174.0344      0.0180
+        Interpolate velo onto osolve   4463.3246    190.1127     16.0783
+                                                                        ov%pressure:     -0.2091E+00 0.0000E+00
+                                                                        osolve%pressure: -0.1797E+00 0.0000E+00
+                      Find bad faces   4493.6690    220.4571     30.3444
+                           Find void   4495.6344    222.4225      1.9654
+                                                                        whole leaf in fluid   196608
+          Define boundary conditions   4495.7223    222.5104      0.0879
+                         wsmp setup1   4495.7450    222.5331      0.0227
+                                                                        n =   699003
+                                                                        nz= 27032364
+                         wsmp setup2   4500.9421    227.7303      5.1972
+ -----------------------------------
+ start of non-linear iteratio      1   4501.1742
+ -----------------------------------
+                        build system   4501.1744      0.0002      0.0002
+                                                                        nelem per proc (min/max)       27063       42687
+                                                                        viscosity range   1.00000E-04  4.24312E+01
+                                                                        viscosity capped in    213008
+                          wsmp solve   4548.5879     47.4137     47.4135
+                                                                        =================================
+                                                                        wsmp: initialisation time:         0.7758930 s
+                                                                        wsmp: ordering time:              22.8445249 s
+                                                                        wsmp: symbolic fact. time:         4.0882258 s
+                                                                        wsmp: Cholesky fact. time:       183.1621981 s
+                                                                        wsmp: Factorisation            14455.4282402 Mflops
+                                                                        wsmp: back substitution time:      1.6904209 s
+                                                                        wsmp: from b to rhs vector:        0.0376720 s
+                                                                        ---------------------------------
+                                                                        wsmp: total solve time:          212.6002140 s
+                                                                        =================================
+                                                                        u :  0.00000E+00  0.10000E+01
+                                                                        v : -0.11412E+00  0.11709E+00
+                                                                        w : -0.15809E+00  0.40457E+00
+                                                                        =================================
+                 Calculate pressures   4761.4766    260.3024    212.8887
+                                                                        raw    pressures : -0.37659E+01  0.36233E+01
+                 Smoothing pressures   4761.8067    260.6325      0.3300
+                                                                        smooth pressures : -0.17019E+00  0.60956E-05
+                do leaf measurements   4761.8964    260.7222      0.0897
+       compute convergence criterion   4762.2554    261.0812      0.3590
+                                                                        velocity_diff_norm = 0.0394506 (tol = 0.0100000)
+                                                                        increase_current_level -> T
+                                                                        current_level -> 4
+                    slicing the cube   4762.2709    261.0967      0.0155
+Compute isostasy and adjust vertical   4762.2715    261.0973      0.0006
+                  Compute divergence   4762.2716    261.0974      0.0001
+                        wsmp_cleanup   4762.5030    261.3288      0.2315
+              Reset surface geometry   4762.5142    261.3400      0.0111
+ -----------------------------------------------------------------------
+           Temperature calculations    4762.5150    261.3408      0.0009
+ -----------------------------------------------------------------------
+                         wsmp setup1   4762.5348    261.3606      0.0197
+                                                                        n =   233001
+                                                                        nz=  5929525
+                         wsmp setup2   4763.5787    262.4045      1.0439
+                        build system   4763.6567    262.4825      0.0779
+                                                                        nelem per proc (min/max)       27058       42687
+                                                                        ===================================
+                                                                        wsmp: initialisation time:         0.0000970 s
+                                                                        wsmp: analysis time:               7.5549400 s
+                                                                        wsmp: nb of exp. non-z in factors   180599703
+                                                                        wsmp: nb of exp. Mflops in factorisation 197655.4935424
+                                                                        wsmp: LU time:                    15.9145639 s
+                                                                        wsmp: actual nb of non-z in factors   180599703
+                                                                        wsmp: actual nb of Mflops in factorisation  12422.9881916
+                                                                        wsmp: back substitution time:      0.2212951 s
+                                                                        wsmp: from b to rhs vector:        0.0123291 s
+                                                                        ---------------------------------
+                                                                        wsmp: total solve time:           24.4434860 s
+                                                                        =================================
+                                                                        temp :  0.00000E+00  0.17149E+01
+                                                                        =================================
+                        wsmp_cleanup   4794.9219    293.7477     31.2653
+                 Advect the surfaces   4795.1561    293.9819      0.2341
+                   Advect surface  1   4795.1562    293.9820      0.0002
+                                                                        S. 1:    1 flips in delaunay2
+                                                                        S. 1:    0 flips in delaunay2
+                    Erode surface  1   4795.2047    294.0305      0.0485
+                 Update cloud fields   4809.2908    308.1166     14.0861
+                    Advect the cloud   4815.3027    314.1285      6.0120
+                        Write output   4850.5740    349.3998     35.2713
+                    End of time step   4862.9956    361.8214     12.4215
+ ----------------------------------------------------------------------- 
+ ----------------------------------------------------------------------- 
+               Start of Step      10   4862.9958
+ ----------------------------------------------------------------------- 
+ ----------------------------------------------------------------------- 
+                  Loop over surfaces   4862.9958      0.0001      0.0001
+                          surface 01   4862.9959      0.0001      0.0000
+                                                                        S. 1:     4356points
+                      refine surface   4862.9959      0.0001      0.0000
+                                                                        S. 1:   0 added ptcls in refine_surface
+                      check delaunay   4863.0113      0.0155      0.0154
+                                                                        S. 1:    0 flips in delaunay2
+                                                                        S. 1:     4356points
+ ----------------------------------------------------------------------- 
+          Start of iteration       1   4863.0145
+ ----------------------------------------------------------------------- 
+                 Create octree solve   4863.0147      0.0001      0.0001
+                                                                        (1) params%griditer < 0                F
+                                                                        (2) current_level==params%levelmax_oct F
+                                                                        (3) increase_current_level             T
+                                                                        converge = (1) & (2) & (3) ->  F
+        compute ref. criterion on ov   4863.0163      0.0018      0.0016
+                                                                        osolve%nleaves=    4096
+        improve osolve for criterion   4863.1603      0.1457      0.1440
+                                                                        max(crit)        =   6.399756E+01
+                                                                        max(crit,w_l_i_f)=   4.596288E+01
+                                                                        osolve%nleaves=    4096
+    Improve osolve for imposed boxes   4863.2831      0.2686      0.1228
+                                                                        osolve%nleaves=  434176
+    Smoothen octree after refinement   4863.5961      0.5815      0.3130
+                                                                        osolve%nleaves=  434176
+                   Build osolve icon   4865.1186      2.1041      1.5226
+             Imbed surface in osolve   4866.8534      3.8389      1.7348
+                embedding surface  1   4866.8538      3.8392      0.0003
+                                                                        osolve%nleaves=  434176
+          Assessing octree stability   5008.0466    145.0321    141.1929
+                                                                        current refine level     4
+                                                                        After criterion based refinement:
+                                                                        nleaves_old=   434176 leaves
+                                                                        nleaves_new=   434176 leaves
+                                                                        After surfaces embedding:
+                                                                        nleaves_old=   434176 leaves
+                                                                        nleaves_new=   434176 leaves
+                                                                        C2: authorise increase of refine level =  T
+                     Change 3D cloud   5008.0602    145.0457      0.0136
+        Interpolate velo onto osolve   5025.8793    162.8648     17.8191
+                                                                        ov%pressure:     -0.1982E+00 0.2447E-04
+                                                                        osolve%pressure: -0.1702E+00 0.6096E-05
+                      Find bad faces   5084.0515    221.0369     58.1721
+                           Find void   5086.0020    222.9875      1.9506
+                                                                        whole leaf in fluid   198134
+          Define boundary conditions   5086.1009    223.0864      0.0989
+                         wsmp setup1   5086.1259    223.1114      0.0250
+                                                                        n =   704745
+                                                                        nz= 27256176
+                         wsmp setup2   5091.3467    228.3322      5.2208
+ -----------------------------------
+ start of non-linear iteratio      1   5091.5777
+ -----------------------------------
+                        build system   5091.5779      0.0002      0.0002
+                                                                        nelem per proc (min/max)       27321       42930
+                                                                        viscosity range   1.00000E-04  2.18906E+02
+                                                                        viscosity capped in    214779
+                          wsmp solve   5139.5141     47.9364     47.9362
+                                                                        =================================
+                                                                        wsmp: initialisation time:         0.7800322 s
+                                                                        wsmp: ordering time:              21.2654810 s
+                                                                        wsmp: symbolic fact. time:         4.0871730 s
+                                                                        wsmp: Cholesky fact. time:       184.8293519 s
+                                                                        wsmp: Factorisation            14603.8460771 Mflops
+                                                                        wsmp: back substitution time:      1.6411312 s
+                                                                        wsmp: from b to rhs vector:        0.0379939 s
+                                                                        ---------------------------------
+                                                                        wsmp: total solve time:          212.6424489 s
+                                                                        =================================
+                                                                        u :  0.00000E+00  0.10000E+01
+                                                                        v : -0.11115E+00  0.12429E+00
+                                                                        w : -0.14192E+00  0.41303E+00
+                                                                        =================================
+                 Calculate pressures   5352.4171    260.8393    212.9029
+                                                                        raw    pressures : -0.34136E+01  0.32712E+01
+                 Smoothing pressures   5352.7484    261.1707      0.3314
+                                                                        smooth pressures : -0.16550E+00  0.53384E-05
+                do leaf measurements   5352.8387    261.2610      0.0903
+       compute convergence criterion   5353.1899    261.6122      0.3512
+                                                                        velocity_diff_norm = 0.0295352 (tol = 0.0100000)
+                                                                        increase_current_level -> T
+                                                                        current_level -> 4
+                    slicing the cube   5353.2056    261.6279      0.0157
+Compute isostasy and adjust vertical   5353.2058    261.6280      0.0002
+                  Compute divergence   5353.2058    261.6281      0.0000
+                        wsmp_cleanup   5353.4312    261.8535      0.2254
+              Reset surface geometry   5353.4436    261.8659      0.0124
+ -----------------------------------------------------------------------
+           Temperature calculations    5353.4445    261.8667      0.0008
+ -----------------------------------------------------------------------
+                         wsmp setup1   5353.4641    261.8864      0.0197
+                                                                        n =   234915
+                                                                        nz=  5978623
+                         wsmp setup2   5354.5182    262.9405      1.0541
+                        build system   5354.6000    263.0223      0.0818
+                                                                        nelem per proc (min/max)       27318       42929
+                                                                        ===================================
+                                                                        wsmp: initialisation time:         0.0000880 s
+                                                                        wsmp: analysis time:               7.4740741 s
+                                                                        wsmp: nb of exp. non-z in factors   182818259
+                                                                        wsmp: nb of exp. Mflops in factorisation 200716.2010137
+                                                                        wsmp: LU time:                    15.9043021 s
+                                                                        wsmp: actual nb of non-z in factors   182818259
+                                                                        wsmp: actual nb of Mflops in factorisation  12623.7221740
+                                                                        wsmp: back substitution time:      0.2345459 s
+                                                                        wsmp: from b to rhs vector:        0.0124640 s
+                                                                        ---------------------------------
+                                                                        wsmp: total solve time:           24.3561862 s
+                                                                        =================================
+                                                                        temp :  0.00000E+00  0.17623E+01
+                                                                        =================================
+                        wsmp_cleanup   5385.8029    294.2252     31.2029
+                 Advect the surfaces   5386.0354    294.4577      0.2325
+                   Advect surface  1   5386.0356    294.4578      0.0002
+                                                                        S. 1:    1 flips in delaunay2
+                                                                        S. 1:    0 flips in delaunay2
+                    Erode surface  1   5386.0842    294.5064      0.0486
+                 Update cloud fields   5400.0519    308.4742     13.9678
+                    Advect the cloud   5406.1550    314.5772      6.1030
+                        Write output   5441.6072    350.0295     35.4523
+                    End of time step   5453.8097    362.2319     12.2024
+ ----------------------------------------------------------------------- 
+ ----------------------------------------------------------------------- 
+               Start of Step      11   5453.8098
+ ----------------------------------------------------------------------- 
+ ----------------------------------------------------------------------- 
+                  Loop over surfaces   5453.8099      0.0001      0.0001
+                          surface 01   5453.8099      0.0001      0.0000
+                                                                        S. 1:     4356points
+                      refine surface   5453.8100      0.0001      0.0000
+                                                                        S. 1:   0 added ptcls in refine_surface
+                      check delaunay   5453.8155      0.0057      0.0056
+                                                                        S. 1:    0 flips in delaunay2
+                                                                        S. 1:     4356points
+ ----------------------------------------------------------------------- 
+          Start of iteration       1   5453.8191
+ ----------------------------------------------------------------------- 
+                 Create octree solve   5453.8192      0.0002      0.0002
+                                                                        (1) params%griditer < 0                F
+                                                                        (2) current_level==params%levelmax_oct F
+                                                                        (3) increase_current_level             T
+                                                                        converge = (1) & (2) & (3) ->  F
+        compute ref. criterion on ov   5453.8208      0.0018      0.0016
+                                                                        osolve%nleaves=    4096
+        improve osolve for criterion   5453.9613      0.1423      0.1405
+                                                                        max(crit)        =   6.399973E+01
+                                                                        max(crit,w_l_i_f)=   4.744129E+01
+                                                                        osolve%nleaves=    4096
+    Improve osolve for imposed boxes   5454.0856      0.2665      0.1243
+                                                                        osolve%nleaves=  434176
+    Smoothen octree after refinement   5454.3987      0.5796      0.3131
+                                                                        osolve%nleaves=  434176
+                   Build osolve icon   5455.9190      2.1000      1.5203
+             Imbed surface in osolve   5457.6357      3.8166      1.7167
+                embedding surface  1   5457.6359      3.8168      0.0002
+                                                                        osolve%nleaves=  434176
+          Assessing octree stability   5669.8114    215.9923    212.1755
+                                                                        current refine level     4
+                                                                        After criterion based refinement:
+                                                                        nleaves_old=   434176 leaves
+                                                                        nleaves_new=   434176 leaves
+                                                                        After surfaces embedding:
+                                                                        nleaves_old=   434176 leaves
+                                                                        nleaves_new=   434176 leaves
+                                                                        C2: authorise increase of refine level =  T
+                     Change 3D cloud   5669.8268    216.0078      0.0154
+        Interpolate velo onto osolve   5686.6203    232.8012     16.7935
+                                                                        ov%pressure:     -0.1927E+00 0.2157E-04
+                                                                        osolve%pressure: -0.1655E+00 0.5338E-05
+                      Find bad faces   5695.8924    242.0734      9.2721
+                           Find void   5697.8359    244.0168      1.9435
+                                                                        whole leaf in fluid   199970
+          Define boundary conditions   5697.9319    244.1128      0.0960
+                         wsmp setup1   5697.9626    244.1435      0.0307
+                                                                        n =   710925
+                                                                        nz= 27503508
+                         wsmp setup2   5703.2375    249.4185      5.2750
+ -----------------------------------
+ start of non-linear iteratio      1   5703.4670
+ -----------------------------------
+                        build system   5703.4671      0.0002      0.0002
+                                                                        nelem per proc (min/max)       27582       43184
+                                                                        viscosity range   1.00000E-04  7.85286E+02
+                                                                        viscosity capped in    216764
+                          wsmp solve   5750.9871     47.5201     47.5200
+                                                                        =================================
+                                                                        wsmp: initialisation time:         0.9279749 s
+                                                                        wsmp: ordering time:              21.5467000 s
+                                                                        wsmp: symbolic fact. time:         4.0976081 s
+                                                                        wsmp: Cholesky fact. time:       197.9714620 s
+                                                                        wsmp: Factorisation            14147.0125909 Mflops
+                                                                        wsmp: back substitution time:     29.5101249 s
+                                                                        wsmp: from b to rhs vector:        0.0382011 s
+                                                                        ---------------------------------
+                                                                        wsmp: total solve time:          254.0933461 s
+                                                                        =================================
+                                                                        u :  0.00000E+00  0.10000E+01
+                                                                        v : -0.14721E+00  0.14765E+00
+                                                                        w : -0.12729E+00  0.43627E+00
+                                                                        =================================
+                 Calculate pressures   6005.3387    301.8717    254.3516
+                                                                        raw    pressures : -0.32791E+01  0.31285E+01
+                 Smoothing pressures   6005.6753    302.2084      0.3367
+                                                                        smooth pressures : -0.16346E+00  0.68000E-06
+                do leaf measurements   6005.7667    302.2997      0.0914
+       compute convergence criterion   6006.1230    302.6560      0.3563
+                                                                        velocity_diff_norm = 0.0229300 (tol = 0.0100000)
+                                                                        increase_current_level -> T
+                                                                        current_level -> 4
+                    slicing the cube   6006.1389    302.6720      0.0159
+Compute isostasy and adjust vertical   6006.1394    302.6724      0.0005
+                  Compute divergence   6006.1394    302.6725      0.0000
+                        wsmp_cleanup   6006.3681    302.9011      0.2287
+              Reset surface geometry   6006.3798    302.9128      0.0117
+ -----------------------------------------------------------------------
+           Temperature calculations    6006.3806    302.9137      0.0009
+ -----------------------------------------------------------------------
+                         wsmp setup1   6006.4006    302.9337      0.0200
+                                                                        n =   236975
+                                                                        nz=  6032899
+                         wsmp setup2   6007.4626    303.9956      1.0620
+                        build system   6007.5443    304.0774      0.0817
+                                                                        nelem per proc (min/max)       27582       43184
+                                                                        ===================================
+                                                                        wsmp: initialisation time:         0.0000839 s
+                                                                        wsmp: analysis time:               7.5835309 s
+                                                                        wsmp: nb of exp. non-z in factors   186032823
+                                                                        wsmp: nb of exp. Mflops in factorisation 207555.3009769
+                                                                        wsmp: LU time:                    16.2542520 s
+                                                                        wsmp: actual nb of non-z in factors   186032823
+                                                                        wsmp: actual nb of Mflops in factorisation  12772.7782992
+                                                                        wsmp: back substitution time:      0.2416210 s
+                                                                        wsmp: from b to rhs vector:        0.0126932 s
+                                                                        ---------------------------------
+                                                                        wsmp: total solve time:           24.9921432 s
+                                                                        =================================
+                                                                        temp :  0.00000E+00  0.18066E+01
+                                                                        =================================
+                        wsmp_cleanup   6039.2890    335.8220     31.7447
+                 Advect the surfaces   6039.5224    336.0554      0.2334
+                   Advect surface  1   6039.5225    336.0556      0.0002
+                                                                        S. 1:    0 flips in delaunay2
+                    Erode surface  1   6039.5700    336.1030      0.0474
+                 Update cloud fields   6053.8811    350.4141     14.3111
+                    Advect the cloud   6060.0835    356.6165      6.2024
+                        Write output   6095.7891    392.3221     35.7056
+                    End of time step   6108.2745    404.8075     12.4854
+ ----------------------------------------------------------------------- 
+ ----------------------------------------------------------------------- 
+               Start of Step      12   6108.2746
+ ----------------------------------------------------------------------- 
+ ----------------------------------------------------------------------- 
+                  Loop over surfaces   6108.2747      0.0001      0.0001
+                          surface 01   6108.2747      0.0001      0.0000
+                                                                        S. 1:     4356points
+                      refine surface   6108.2748      0.0001      0.0000
+                                                                        S. 1:   0 added ptcls in refine_surface
+                      check delaunay   6108.2810      0.0063      0.0062
+                                                                        S. 1:    0 flips in delaunay2
+                                                                        S. 1:     4356points
+ ----------------------------------------------------------------------- 
+          Start of iteration       1   6108.2842
+ ----------------------------------------------------------------------- 
+                 Create octree solve   6108.2844      0.0002      0.0002
+                                                                        (1) params%griditer < 0                F
+                                                                        (2) current_level==params%levelmax_oct F
+                                                                        (3) increase_current_level             T
+                                                                        converge = (1) & (2) & (3) ->  F
+        compute ref. criterion on ov   6108.2860      0.0018      0.0016
+                                                                        osolve%nleaves=    4096
+        improve osolve for criterion   6108.4262      0.1420      0.1402
+                                                                        max(crit)        =   6.399997E+01
+                                                                        max(crit,w_l_i_f)=   4.871429E+01
+                                                                        osolve%nleaves=    4096
+    Improve osolve for imposed boxes   6108.5503      0.2661      0.1241
+                                                                        osolve%nleaves=  434176
+    Smoothen octree after refinement   6108.8666      0.5824      0.3163
+                                                                        osolve%nleaves=  434176
+                   Build osolve icon   6110.4121      2.1278      1.5454
+             Imbed surface in osolve   6112.1853      3.9010      1.7732
+                embedding surface  1   6112.1854      3.9012      0.0002
+                                                                        osolve%nleaves=  434176
+          Assessing octree stability   6266.0075    157.7232    153.8220
+                                                                        current refine level     4
+                                                                        After criterion based refinement:
+                                                                        nleaves_old=   434176 leaves
+                                                                        nleaves_new=   434176 leaves
+                                                                        After surfaces embedding:
+                                                                        nleaves_old=   434176 leaves
+                                                                        nleaves_new=   434176 leaves
+                                                                        C2: authorise increase of refine level =  T
+                     Change 3D cloud   6266.0223    157.7381      0.0149
+        Interpolate velo onto osolve   6285.5921    177.3078     19.5697
+                                                                        ov%pressure:     -0.1899E+00 0.4187E-05
+                                                                        osolve%pressure: -0.1635E+00 0.6800E-06
+                      Find bad faces   6328.0884    219.8042     42.4964
+                           Find void   6330.0631    221.7788      1.9747
+                                                                        whole leaf in fluid   200686
+          Define boundary conditions   6330.1640    221.8798      0.1009
+                         wsmp setup1   6330.1937    221.9095      0.0297
+                                                                        n =   713055
+                                                                        nz= 27590838
+                         wsmp setup2   6335.5008    227.2165      5.3070
+ -----------------------------------
+ start of non-linear iteratio      1   6335.7356
+ -----------------------------------
+                        build system   6335.7358      0.0002      0.0002
+                                                                        nelem per proc (min/max)       27672       43259
+                                                                        viscosity range   1.00000E-04  1.00000E+03
+                                                                        viscosity capped in    217474
+                          wsmp solve   6384.3182     48.5826     48.5824
diff --git a/test/input.txt b/test/input.txt
new file mode 100644
index 0000000000000000000000000000000000000000..7c042c9736e9052d62773d61d3029e3730be046f
--- /dev/null
+++ b/test/input.txt
@@ -0,0 +1,625 @@
+!------------------------------------------------------------------------------|
+!------------------------------------------------------------------------------|
+!                                                                              |
+!              ||===\\                                                         | 
+!              ||    \\                                                        |
+!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
+!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
+!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
+!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
+!                                                                              |
+!              Input File                                                      |
+!------------------------------------------------------------------------------|
+!------------------------------------------------------------------------------|
+
+This file is read by two subroutines in the code:
+- read_controlling_parameters
+- read_input_file
+The indications between [] indicate the type of the read parameter. It can be
+an integer [int], a real*8 [dp], a character chain [char] or a boolean [bool].
+
+
+CONTROLLING PARAMETERS
+_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
+
+[int]debug is a parameter that switches on/off various prints and outputs
+(the level of printing for error, warning and messages of the solver id%cntl(4)
+is set to the debug value).
+if debug is zero, no debugging 
+if debug is equal to one, this triggers the terminal display of some key parameters 
+if debug is equal to two, same as debug equal to one, and the code produces various 
+vtk files in the DEBUG subdirectories (surfaces, olsf, ...), as well as cross sections 
+for all nonlinear iterations. Careful, memory consuming !
+
+
+      debug = 1
+
+[bool]doDoRuRe is a flag that triggers the production of output files needed to produce
+the DoRuRes. DoRuRe stands for 'Douar Run Report'. 
+
+       doDoRuRe = F
+
+[bool]compute_qp_gram triggers the production of qpgrams for every grid.
+
+      compute_qpgram = F
+
+[bool]compute_reaction_forces toggles on/off the reaction forces computations.
+
+      compute_reaction_forces = F
+
+RESTART
+_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
+
+[int]irestart is a restart flag; if irestart is not 0, the run will restart from
+an output file given by [char]restartfile and at step.
+
+      irestart = 0
+
+      restartfile=OUT/time_0187.bin
+
+
+TIMESTEPPING
+_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
+
+[dp]dt is the time step length (if dt is negative, courant conidition is used and
+automatic time stepping is turned on)
+
+      dt=1.d-3
+
+[int]nstep is the number of time steps
+
+      nstep = 200
+
+[dp]courant is only used when dt is negative; it determines the size of the time
+step from the maximum value of the velocity field amplitude. The time step s
+the product of courant by the ratio of the smallest leaf size by the maximum
+velocity.
+
+      courant=.5d0
+
+[bool]normaladvect is a flag used to determine which algorithm to use to calculate
+the new geoletry of the normals to the surfaces at the nodes on the surfaces
+if normaladvect is T, the normals are advected using the velocity gradient
+if normaladvect is F, the normals are re-computed from the geometry of the
+surface
+
+      normaladvect = T
+
+GRID ITERATIONS
+_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
+
+[int]griditer is a flag that allows for nonlinear iterations; when positive, a fixed number
+(griditer) of iterations is permitted; when negative, the number of
+nonlinear iterations is determined by a convergence criterion.
+
+      griditer = 1
+
+[dp]octree_refine_ratio is the threshold value used to determine whether the octree 
+has converged or not. the larger the value, the less stringent the test.
+
+      octree_refine_ratio=1.025d0
+
+
+NONLINEAR ITERATIONS
+_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
+
+
+[int]nonlinear_iterations is the maximum number of nonlinear iterations (i.e. the 
+iterations on a given constant grid)
+if nonlinear_iterations is positive, it simply is the number of nonlinear
+iterations performed for each grid. When negative it indicates an upper bound
+of nonlinear iterations, but the actual number of nonlinear iterations is
+determined by a convergence criterion (see the 'tol' parameter)
+
+      nonlinear_iterations = 1
+
+[dp]tol is the relative tolerance used to estimate convergence on the computed
+velocity field
+
+      tol=0.01d0
+
+[bool]adaptive_tol is a flag that toggles on/off the evolution of the tol parameter
+with the grid level: when velocity convergence is reached on a grid, the
+latest meaure of the velocity difference between the two last obtained
+solutions is put in tol, so that on the following generated grid, the solution
+reaches at least the same level of convergence. It allows to start with a not
+too stringent value of tol at uniform octree level that evolves with every
+grid, assuming that increasing the level of refinement of the octree allows to
+better capture the solution, hence allowing a tighter convergence. 
+
+      adaptive_tol = F
+
+OCTREES
+_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
+
+[int]leveluniform_oct is the level of uniform discretization of space; note that a
+level is a power of two used to divide the unit cube
+
+      leveluniform_oct = 4
+
+[int]levelmax_oct is maximum level of octree discretization
+
+      levelmax_oct = 6
+
+[bool]ismooth is a flag to impose an additional level of smoothing after refinement
+for the surfaces and strain rate. It ensures that no leaf is flanked by
+other leaves differing by more than 1 level of refinement
+
+      ismooth = F
+
+[int]noctreemax is the maximum size of any octree used in all computations
+
+      noctreemax=10000000
+
+[dp]refine_ratio is used to determine octree refinement based on a given criterion. 
+All leaves where the criterion is larger than refine_ratio times the maximum 
+of this criterion are refined 
+
+      refine_ratio=-10000
+
+[int]refine_criterion determines which refinement algorithm is to be used.
+Several criteria exist for the refinement of the osolve octree. 1 is the second 
+invariant of the deviatoric strain-rate tensor; 2 is the sum of the squares of the
+diagonal terms of the deviatoric strain-rate tensor; 3 is the second invariant
+of the deviatoric strain rate tensor timses the leaf size. any other value sets
+the criterion to zero and leads to no refinement. 
+
+      refine_criterion=1
+
+[int]initial_refine_level is the initial level at which the refinement of the
+octree will be performed. it has to be smaller than levelmax_oct
+this should be used (different from levelmax_oct) in case the flow is very 
+localized (nonlinear/plastic analysis)
+
+      initial_refine_level = 4
+
+[bool]renumber_nodes is a flag that can toggle on/off the renumbering of the nodes
+by mean of Sloan's algorithm (T/F)
+
+      renumber_nodes = T
+
+PRESSURE
+_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
+[int]smoothing_type is a parameter allows to choose which type of smoothing is to 
+be applied to the pressure field: 0 is none, 1 is center->nodes->center, 2 is the 
+same, but weighted by neighbouring elemental volumes, 3 is regular grid+SPH, and 4 is SPH.  
+
+      smoothing_type = 1
+
+CLOUD
+_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
+
+[int]npmin and [int]npmax are used to update the 3D volumetric cloud. npmin corresponds
+to the minimum number of particles in any leaf; npmax is the maximum allowable
+number in any leaf
+
+      npmin=8
+      npmax=27
+
+
+FEM + DIVFEM + MUMPS
+_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
+
+[int]levelcut is the number of levels used to estimate the volume integrals in the
+divFEM algorithmi for cut cells; by testing, it has been estimated that a value
+of 2 is very accurate
+
+      levelcut=2
+
+[int]levelapprox is the number of levels used to estimate the remaining integrals
+using an improved version of Marthijn's clever algorithm, usually 3 is plenty
+
+      levelapprox=3
+
+[dp]penalty is a global penalty parameter used to impose the bad faces or
+incompatible faces linear constraints
+
+      penalty=1.d8
+
+[bool]excl_vol is a parameter that toggles off the assumption that lsf's are built
+on top of one another
+
+      excl_vol = F 
+
+
+TEMPERATURE
+_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
+[dp] ztemp is the height interval between which a linear temperature gradient 
+is set: the temperature is 1 at the bottom, and 0 at ztemp.
+[dp] tempscale
+
+      calculate_temp = T
+
+      ztemp = .1d0
+
+      tempscale=250.d0
+
+MATERIALS
+_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
+
+[int]nmat is number of materials
+
+      nmat=1
+
+[int]material0 determines what is the material above the 1st surface (the free
+surface) if material0 is 0 then it is the void and the properties "0" are used
+for that part of the model; otherwise the material is one of the material,
+comprised between 1 and nmat
+
+      material0 = 0
+
+[dp]densityi, [dp]viscosityi and [dp]penaltyi are the density, viscosity and
+incompressibility used for material i; there should be nmat sets of material
+properties; there should also be a nil material if material0 has been set to 0
+[dp]expon is the nonlinear viscosity exponent
+[dp]diffusivity is the heat diffusivity
+[dp]heat is the heat production
+
+[char]plasticity_type is the type of plasticity
+- No    -> no plasticity, purely (nonlinear) viscous material
+- vM    -> von Mises yield criterion
+        -> [dp]plasticity_1st_param is the yield value
+- DPI   -> Drucker-Prager type of yield criterion
+        -> the yield locus passes through the outer apices of the Mohr-Coulomb hexagon
+        -> [dp]plasticity_1st_param is the angle phi
+        -> [dp]plasticity_2nd_param is the cohesion c
+- DPII  -> Drucker-Prager type of yield criterion
+        -> the yield locus passes through the inner apices of the Mohr-Coulomb hexagon
+        -> [dp]plasticity_1st_param is the angle phi
+        -> [dp]plasticity_2nd_param is the cohesion c
+- DPIII -> Drucker-Prager type of yield criterion
+        -> [dp]plasticity_1st_param is alpha
+        -> [dp]plasticity_2nd_param is k
+- MC    -> Mohr-Coulomb type of yield criterion
+        -> [dp]plasticity_1st_param is the angle phi
+        -> [dp]plasticity_2nd_param is the cohesion c
+
+      density0              = 0.d0
+      viscosity0            = 1.d-5
+      penalty0              = 1.d8
+      expon0                = 1.d0
+      diffusivity0          = 1.d0
+      heat0                 = 0.d0
+      activationenergy0     = 0.d0
+      plasticity_type0      = MC
+      plasticity_1st_param0 = 15.d0
+      plasticity_2nd_param0 = 1.72d-3
+      plasticity_3rd_param0 = 0.5d0
+      plasticity_4th_param0 = 1.5d0
+      plasticity_5th_param0 = 5.d0
+
+      density1              = -.83d0
+      viscosity1            = 9.136d-8
+      penalty1              = 1.d8
+      expon1                = 4.d0
+      diffusivity1          = 1.25d-2
+      heat1                 = 0.d0
+      activationenergy1     = 223.d3
+      plasticity_type1      = MC
+      plasticity_1st_param1 = 15.d0
+      plasticity_2nd_param1 = 10.d-3
+old value for plasticity_2nd_param1 is 3.13d-3
+      plasticity_3rd_param1 = 0.5d0
+      plasticity_4th_param1 = 1.5d0
+      plasticity_5th_param1 = 5.d0
+
+[dp]viscositymin and viscositymax are bounds on the viscosity (if negative bound is not imposed)
+These bounds are introduced to prevent the viscosity to reach unrealistic values, especially when
+using non-linear (power-law or brittle) rheologies
+
+      viscositymin=1.d-4
+      viscositymax=1.d3
+
+SURFACES
+_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
+
+[int]ns is number of surfaces to track
+
+      ns=1
+
+for each surface, one needs to define a levelt, itype, material, surface_type, activation_time 
+and surface_params .
+- [int]levelt is the inital level for the particles on the surface; to be accurate
+and avoid wholes in the surface during definition of the lsf, one should use
+levelt eq to levelmax_oct+1 for all surfaces as a minimum value; 
+- [int]itype should be 1 for foldable surfaces or 0 for nonfoldable surfaces; 
+- [int]material is the material type refering to the table of material available (max nmat); 
+- [dp]activation_time is the time the surface becomes active (before that time, it
+  is glued to the 0th surface). this parameter is useful when defining stratigraphic horizons; 
+  default is -1, ie the surface is not glued to the free surface
+- [int]surface_type is comprised between 1 and 8. 
+  1 corresponds to  a flat surface, 
+    -> surface_param_01 is the z level
+  2 to a rectangular emboss, 
+    -> surface_param_01 is the z level
+    -> surface_param_02 and 03 are x1,x2 
+    -> surface_param_04 and 05 are y1,y2 
+    -> surface_param_06 is the thickness 
+  3 to a convex spherical emboss, 
+    -> surface_param_01 is the z level
+    -> surface_param_02 and 03 are x0,y0 
+    -> surface_param_04 is the radius 
+  4 to concave spherical emboss,
+    -> surface_param_01 is the z level
+    -> surface_param_02 and 03 are x0,y0 
+    -> surface_param_04 is the radius 
+  5 to a double rectangular emboss, 
+    -> surface_param_01 is the z level
+    -> surface_param_02 and 03 are x1,x2 
+    -> surface_param_04 and 05 are x3,x4 
+    -> surface_param_06 and 07 are y1,y2 
+    -> surface_param_08 and 09 are y3,y4 
+    -> surface_param_10 is the thickness 
+  6 to a sinus, 
+    -> surface_param_01 is the z level
+    -> surface_param_02 is the wavelength 
+    -> surface_param_03 is the amplitude 
+  7 to a noisy surface, 
+    -> surface_param_01 is the z level
+    -> surface_param_02 is the noise amplitude
+  8 to a double sinus.
+    -> surface_param_01 is the z level
+    -> surface_param_02 is the x-wavelength 
+    -> surface_param_03 is the x-amplitude 
+    -> surface_param_04 is the y-wavelength 
+    -> surface_param_05 is the y-amplitude 
+- [int]leveloct is the level at which the octree will be refined in the vicinity of the surface.
+
+      levelt1             = 6
+      itype1              = 0
+      surface_type_1      = 1
+      rand1               = T
+      surface_param_01_1  = .1d0
+      material1           = 1
+      activation_time_1   = -1. 
+      leveloct1           = 6
+      stretch1            = 1.5d0
+      anglemax1           = 180.d0
+      criterion1          = 1
+      anglemaxoctree1     = 180.d0
+      spread_surface_points1 = 1
+
+[int]niter_move is the number of iterations used to update particle positions
+using an implicit, mid-point algorithm (default is 10)
+
+      niter_move = 10
+
+[dp]stretch is the maximum allowed increase in linear length between two initially
+adjacent particles on any surface; when this stretch is achieved, a new
+particle is inserted on the surface, half-way along the stretched edge
+
+      stretch = 1.5d0
+
+[dp]anglemax is the maximum allowed angle between two adjacent normals
+when the angle is reached a new point is inserted beteen the two points to
+reduce the angle between the two normals
+
+      anglemax=180.d0
+
+
+[int]criterion is criterion used to define the octree in the vicinity of the
+sufaces; criterion 1 corresponds to imposing that all leaves that are cut by any of the
+surfaces must be at level levelmax_oct; criterion 2 corresponds
+to imposing that discretization is proportional to the curvature of the
+surface; curvature is calculated from the local divergence of the normals.
+criterion 3 corresponds to imposing that all leaves that contain at
+least one particle of any surface is at levelmax_oct; 
+
+      criterion = 2
+
+[dp]anglemaxoctree is only defined for criterion 2; t is the maximum allowable angle
+between two adjacent normals; if the angle is greater than anglemaxoctree, the local
+octree leaves are forced to be at level levelmax_oct; otherwise they are
+proportionally larger (smaller levels) (default is 10)
+
+     anglemaxoctree = 10.d0
+
+
+REFINEMENT IN BOXES
+_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
+
+[int]nboxes is the number of boxes in which the user imposes a set level of
+discretization
+
+      nboxes = 1
+
+for each box we need to specify the two end corners and the level
+the syntax is [dp]boxnx0, [dp]boxnx1, [dp]boxny0, [dp]boxny1, [dp]boxnz0, [dp]boxnz1
+and boxnlevel where n is the box number
+      box1x0=0.d0
+      box1x1=1.d0
+      box1y0=0.d0
+      box1y1=1.d0
+      box1z0=0.0d0
+      box1z1=0.2d0
+      box1level=7
+
+      box2x0=0.1d0
+      box2x1=0.4d0
+      box2y0=0.0d0
+      box2y1=0.5d0
+      box2z0=0.0d0
+      box2z1=0.2d0
+      box2level=7
+
+      box3x0=0.6d0
+      box3x1=0.9d0
+      box3y0=0.5d0
+      box3y1=1.0d0
+      box3z0=0.0d0
+      box3z1=0.2d0
+      box3level=7
+
+REFINEMENT ON CUBE FACES
+_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
+[bool]ref_on_faces toggles on/off the user imposed refinement on faces.
+for each of the six faces of the cube, on defines the level at which the
+desired area is to be refined. this area is given by bottom, top, left and
+right coordinates that are comprised between 0 and 1
+faces 1,2,3,4,5,6 respectively correspond to x=0,x=1,y=0,y=1,z=0,z=1
+
+      ref_on_faces = F
+
+      level_face1=5
+      b1=.02
+      t1=.51
+      l1=.11
+      r1=.81
+      level_face2=5
+      b2=.45
+      t2=.55
+      l2=.46
+      r2=.56
+      level_face3=5
+      b3=.0
+      t3=.2
+      l3=.3
+      r3=.7
+      level_face4=5
+      b4=.4
+      t4=.5
+      l4=.14
+      r4=.4
+      level_face5=5
+      b5=0.
+      t5=1.
+      l5=0.
+      r5=0.75
+      level_face6=5
+      b6=.26
+      t6=.56
+      l6=.16
+      r6=.86
+
+
+EROSION
+_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
+
+[bool]erosion is a flag that toggles on/off the erosion. 
+
+      erosion = T
+
+if erosion is on, one also needs to set the erosion level/height, [dp]zerosion
+Note that this is a first attempt at erosion; in future versions, DOUAR
+should be easily linked to a surface processes model like CASCADE
+
+      zerosion=.1d0
+
+If erosion is on, one also needs to define a length scale and a velocity scale to
+properly translate the uplift rate produced by DOUAR into something that is adequate
+for CASCADE; 
+[dp]length_scale is the scale of the unit model in km 
+[dp]velocity_scale is the velocity scale in km/Myr. 
+if length_scale is negative, erosion is assumed to be perfect (no call to CASCADE is needed)
+
+      length_scale=200.d0
+      velocity_scale=10.d0
+
+One then needs to define the erosion constants in CASCADE
+[dp]fluvial_erosion is the fluvial erosion constant in 1/m^2
+[dp]diffusion_erosion is the diffusion erosion constant in m^2/yr
+(4d-2 32d-2)
+
+      fluvial_erosion=2.d-2
+      diffusion_erosion=16.d-2
+
+One finally needs to specify the boundary conditions for CASCADE
+if [int]baselevelx0 is set to 1 then the boundary at x=0 is set at baselevel (water and sediment exit)
+if [int]baselevelx1 is set to 1 then the boundary at x=1 is set at baselevel (water and sediment exit)
+if [int]baselevely0 is set to 1 then the boundary at y=0 is set at baselevel (water and sediment exit)
+if [int]baselevely1 is set to 1 then the boundary at y=1 is set at baselevel (water and sediment exit)
+
+      baselevelx0=1
+      baselevelx1=1
+      baselevely0=0
+      baselevely1=0
+
+MATRIX VISUALISATION
+_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
+
+[bool]visualise_matrix allows the user to turn on the visual representation of the 
+matrices used in the code. Be careful, the generated postscript files are huge !
+
+      visualise_matrix = F
+
+CROSS SECTIONS
+_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
+
+[int]nsections is the number of cross-sections to be output. [dp]scale is the scale used to produce the postscript file.
+[int]xyz takes values 1, 2 or 3, and corresponds to planes defined by x=constant, y=constant, and z=constant resp.
+The [bool]flags are self explanatory. [char]colormap is the chosen colormap (jet or hot). [int]ncolours is the number of
+colours used to produce the plot. 
+
+nsections = 0
+
+xyz_1          = 1 
+slice_1        = 0.5001
+flag_press_1   = T
+flag_spress_1  = T 
+flag_e2d_1     = T
+flag_e3d_1     = F
+flag_lode_1    = F 
+flag_crit_1    = F 
+flag_grid_1    = T 
+flag_mu_1      = F   
+flag_u_1       = F
+flag_v_1       = F
+flag_w_1       = F
+flag_q_1       = F 
+flag_uvw_1     = F 
+flag_lsf_1     = F 
+flag_vfield_1  = F 
+flag_colour_1  = T
+flag_plastic_1 = F
+flag_velvect_1 = T
+scale_1        = 500.
+colormap_1     = jet
+ncolours_1     = 256
+
+xyz_2          = 2
+slice_2        = 0.901
+flag_press_2   = T
+flag_e2d_2     = T
+flag_e3d_2     = F 
+flag_lode_2    = F 
+flag_crit_2    = F 
+flag_grid_2    = F 
+flag_mu_2      = F   
+flag_u_2       = F
+flag_v_2       = F
+flag_w_2       = F
+flag_q_2       = F 
+flag_uvw_2     = F 
+flag_lsf_2     = F 
+flag_vfield_2  = F 
+flag_colour_2  = T
+flag_plastic_2 = F
+flag_velvect_2 = F
+scale_2        = 800.
+colormap_2     = jet 
+ncolours_2     = 256
+
+
+xyz_3          = 3
+slice_3        = 0.0010
+flag_press_3   = F 
+flag_e2d_3     = F
+flag_e3d_3     = F
+flag_lode_3    = F 
+flag_crit_3    = F 
+flag_grid_3    = F 
+flag_mu_3      = F   
+flag_u_3       = F
+flag_v_3       = F
+flag_w_3       = F
+flag_q_3       = F
+flag_uvw_3     = F 
+flag_lsf_3     = F 
+flag_vfield_3  = F 
+flag_colour_3  = T
+flag_plastic_3 = F
+flag_velvect_3 = F
+scale_3        = 400.
+colormap_3     = jet
+ncolours_3     = 256
+