Skip to content
Snippets Groups Projects
input.090710.txt 22.8 KiB
Newer Older
  • Learn to ignore specific revisions
  • Douglas Guptill's avatar
    Douglas Guptill committed
    1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718
    !------------------------------------------------------------------------------|
    !------------------------------------------------------------------------------|
    !                                                                              |
    !              ||===\\                                                         |
    !              ||    \\                                                        |
    !              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
    !              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
    !              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
    !              ||===//     \\==//    \\==\\  \\==\\  ||                        |
    !                                                                              |
    !              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_0170.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-4
    
    [int]nstep is the number of time steps
    
         nstep = 500
    
    [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 = -10
    
    [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=.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 = -10
    
    [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 = 5
    
    [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=100000
    
    [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=-5500
    
    [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=0
    
    [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 = 6
    
    [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 = 0
    
    
    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=1
         npmax=4
    
    
    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 = F
    
         ztemp = .357d0
    
         tempscale=1350.d0
    
    
    ISOSTASY
    _-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
    
         isostasy = T
    
         flexure = T
    
    
    MATERIALS
    _-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
    
    [int]nmat is number of materials
    
         nmat=3
    
    [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      = No
    
         density1              = -.83d0
         viscosity1            = 1.d0
         penalty1              = 1.d8
         expon1                = 1.d0
         diffusivity1          = 0.59d0
         heat1                 = 0.d0
         activationenergy1     = 0.d0
         plasticity_type1      = MC
         plasticity_1st_param1 = 15.d0
         plasticity_2nd_param1 = 1.72d-3
         plasticity_3rd_param1 = 0.5d0
         plasticity_4th_param1 = 1.5d0
         plasticity_5th_param1 = 5.d0
    
         density2              = -.83d0
         viscosity2            = 1.d-3
         penalty2              = 1.d8
         expon2                = 1.d0
         diffusivity2          = 0.59d0
         heat2                 = 0.d0
         activationenergy2     = 0.d0
         plasticity_type2      = MC
         plasticity_1st_param2 = 15.d0
         plasticity_2nd_param2 = 1.72d-3
         plasticity_3rd_param2 = 0.5d0
         plasticity_4th_param2 = 1.5d0
         plasticity_5th_param2 = 5.d0
    
         density3              = -1.d0
         viscosity3            = 1.d1
         penalty3              = 1.d8
         expon3                = 1.d0
         diffusivity3          = 0.59d0
         heat3                 = 0.d0
         activationenergy3     = 0.d0
         plasticity_type3      = MC
         plasticity_1st_param3 = 15.d0
         plasticity_2nd_param3 = 1.72d-3 
         plasticity_3rd_param3 = 0.5d0
         plasticity_4th_param3 = 1.5d0
         plasticity_5th_param3 = 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-5
         viscositymax=1.d5
    
    
    SURFACES
    _-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
    
    [int]ns is number of surfaces to track
    
         ns=3
    
    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 holes 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             = 7
         itype1              = 0
         surface_type_1      = 1
         rand1               = T
         surface_param_01_1  = 0.357d0
         material1           = 1
         activation_time_1   = -1. 
         leveloct1           = 6
         stretch1            = 1.5d0
         anglemax1           = 180.d0
         criterion1          = 1
         anglemaxoctree1     = 180.d0
         spread_surface_points1 = 1
    
         levelt2             = 7
         itype2              = 0
         surface_type_2      = 1
         rand2               = T
         surface_param_01_2  = 0.32575d0
         material2           = 2
         activation_time_2   = -1.
         leveloct2           = 6
         stretch2            = 1.5d0
         anglemax2           = 180.d0
         criterion2          = 1
         anglemaxoctree2     = 180.d0
         spread_surface_points2 = 0
    
         levelt3             = 7
         itype3              = 0
         surface_type_3      = 1
         rand3               = T
         surface_param_01_3  = 0.2945d0
         material3           = 3
         activation_time_3   = -1.
         leveloct3           = 6
         stretch3            = 1.5d0
         anglemax3           = 180.d0
         criterion3          = 1
         anglemaxoctree3     = 180.d0
         spread_surface_points3 = 0
    
    [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 = 2
    
    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=.25d0
         box1x1=.75d0
         box1y0=.0d0
         box1y1=.6d0
         box1z0=.2d0
         box1z1=.4d0
         box1level=6
    
         box2x0=0.d0
         box2x1=.75d0
         box2y0=.4d0 
         box2y1=.6d0
         box2z0=.2d0
         box2z1=.4d0
         box2level=6
    
         box3x0=0.d0
         box3x1=1.d0
         box3y0=0.d0
         box3y1=1.d0
         box3z0=0.d0
         box3z1=.001d0
         box3level=6
    
         box4x0=.4875d0
         box4x1=.5125d0
         box4y0=.4875d0
         box4y1=.5125d0
         box4z0=.0d0
         box4z1=1.d0
         box4level=9
    
    
    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=6
         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 = F
    
    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=.357d0
    
    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=350.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=4.d-2
         diffusion_erosion=32.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=1
    
    
    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