diff --git a/src/Makefile.inc.geodyncomp b/src/Makefile.inc.geodyncomp
new file mode 100644
index 0000000000000000000000000000000000000000..66ffe167e517f291c1300dac1c782599ced22c53
--- /dev/null
+++ b/src/Makefile.inc.geodyncomp
@@ -0,0 +1,75 @@
+# Define our compilers.................................
+#
+# on geodyncomp, the defaults are OpenMPI (1.4.3),
+# built with the GNU compilers
+F90 = mpif90
+F77 = mpif90
+cc = mpicc
+
+# we take the standard ar
+AR=ar
+
+# There are no BITS
+BITS=
+
+# no MKL
+# IMKL = /scinet/gpc/intel/Compiler/11.1/056/mkl/include
+# IMKL = /scinet/gpc/intel/Compiler/11.1/072/mkl/include
+# INCLUDE=-I$(IMKL)
+
+# Define FORTRAN compiler flags.......................
+#
+F90FLAGS=-ffree-line-length-0
+
+# Define C compiler flags.............................
+#
+CFLAGS=
+
+# no Intel MKL library
+#
+# LMKL = /scinet/gpc/intel/Compiler/11.1/056/mkl/lib
+# LMKL = /scinet/gpc/intel/Compiler/11.1/072/mkl/lib
+
+# MPI libraries
+# (use the default)
+
+# WSMP Library
+#
+# WSMP = /home/dwhipp/software/wsmp/wsmp-Linux64/lib/Intel
+# WSMP = /home/dwhipp/software/wsmp/wsmp-Linux64-MPI2
+# WSMP = /home/dwhipp/software/wsmp/wsmp-Linux64-IMPI
+WSMP = -L/opt/wsmp/wsmp-Linux64/lib/GNU/openmpi -lpwsmp64
+
+# Not using Google's TCMalloc
+# TCML = /home/dwhipp/software/tcml
+
+LIBS = \
+ -LOCTREE -lOctree$(BITS) \
+ -LNN -lnn_f$(BITS) \
+ -LNN -lnn_c$(BITS) \
+ -LCASCADE -lcascade$(BITS) \
+ -LRESAMPLE -lresample$(BITS) \
+ $(WSMP) -lpthread
+
+# The Makefiles look at compile output with this
+PAGER=more
+
+# compile rules.......................................
+#
+COMPILE_OUT=$(NAME)$(BITS).compile
+.SUFFIXES:
+.SUFFIXES: .o .c .f .f90
+.f90.o:
+ @echo "--" >>$(COMPILE_OUT) 2>&1
+ @echo "$(F90) $(F90FLAGS) $(INCLUDES) -c $<" >>$(COMPILE_OUT) 2>&1
+ $(F90) $(F90FLAGS) $(INCLUDES) -c $< >>$(COMPILE_OUT) 2>&1
+
+.f.o:
+ @echo "--" >>$(COMPILE_OUT) 2>&1
+ @echo "$(F77) $(F77FLAGS) $(INCLUDES) -c $<" >>$(COMPILE_OUT) 2>&1
+ $(F77) $(F77FLAGS) $(INCLUDES) -c $< >>$(COMPILE_OUT) 2>&1
+
+.c.o:
+ @echo "--" >>$(COMPILE_OUT) 2>&1
+ @echo "$(CC) $(CFLAGS) -c $<" >>$(COMPILE_OUT) 2>&1
+ $(CC) $(CFLAGS) -c $< >>$(COMPILE_OUT) 2>&1
diff --git a/src/define_bc_2011-01-21.f90 b/src/define_bc_2011-01-21.f90
new file mode 100644
index 0000000000000000000000000000000000000000..0600ddd9a714c3578174dea3084aa51d31957b08
--- /dev/null
+++ b/src/define_bc_2011-01-21.f90
@@ -0,0 +1,182 @@
+!------------------------------------------------------------------------------|
+!------------------------------------------------------------------------------|
+! |
+! ||===\\ |
+! || \\ |
+! || || //==\\ || || //==|| ||/==\\ |
+! || || || || || || || || || || |
+! || // || || || || || || || |
+! ||===// \\==// \\==\\ \\==\\ || |
+! |
+!------------------------------------------------------------------------------|
+!------------------------------------------------------------------------------|
+! |
+! DEFINE_BC Feb. 2009 |
+! |
+!------------------------------------------------------------------------------|
+!------------------------------------------------------------------------------|
+
+subroutine define_bc (params,osolve,vo,zi)
+
+!------------------------------------------------------------------------------|
+!(((((((((((((((( Purpose of the routine )))))))))))))))))))))))))))))))))))))
+!------------------------------------------------------------------------------|
+! this routine assigns the boundary condition for the Tibet experiment
+
+!------------------------------------------------------------------------------|
+!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
+!------------------------------------------------------------------------------|
+
+use definitions
+
+implicit none
+
+type (parameters) params
+type (octreesolve) osolve
+type (void) vo
+type (ziso) zi
+
+!------------------------------------------------------------------------------|
+!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
+!------------------------------------------------------------------------------|
+
+integer i,iproc,nproc,ierr
+double precision eps,lsf0,pi,lorig,h,x1,x2,phi,spu,spv,yend,cper,cscl,xstart
+double precision ystart,theta,l,vin,vzfluxscl,cntvel,dxy,xend,xsym,ymax,xwidth
+double precision ywidth
+double precision,dimension(:),allocatable :: x0,ldisp
+integer ie,ij,j,jp,nelemx,nelemz
+logical firstcall
+
+!------------------------------------------------------------------------------|
+!------------------------------------------------------------------------------|
+
+INCLUDE 'mpif.h'
+
+call mpi_comm_size (mpi_comm_world,nproc,ierr)
+call mpi_comm_rank (mpi_comm_world,iproc,ierr)
+
+eps=1.d-10
+osolve%kfix=0
+osolve%kfixt=0
+pi=atan(1.d0)*4.d0
+
+select case(trim(params%infile))
+
+case('input.txt','input.small.txt')
+ do i=1,osolve%nnode
+ if (osolve%x(i).lt.eps) then
+ osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=0.d0
+ osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
+ osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
+ endif
+
+ if (osolve%x(i).gt.1.d0-eps) then
+ osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=0.d0
+ osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
+ osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
+ endif
+
+ if (osolve%y(i).lt.eps) then
+! osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=0.d0
+ osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
+! osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
+ endif
+
+ if (osolve%y(i).gt.1.d0-eps) then
+! osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=0.d0
+ osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
+! osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
+ endif
+
+ if (osolve%z(i).lt.eps) then
+ osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=0.d0
+ osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
+ osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
+ osolve%kfixt(i)=1
+ osolve%temp(i)=1.d0
+ endif
+
+ if (osolve%z(i).gt.1.d0-eps) then
+ osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=0.d0
+ osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
+ osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
+ osolve%kfixt(i)=1
+ osolve%temp(i)=0.d0
+ endif
+
+ if (.not.vo%influid(i)) then
+ osolve%kfixt(i)=1
+ osolve%temp(i)=0.d0
+ endif
+
+ ! Set velocity B/C everywhere to disallow flow in the y-direction
+ osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
+ enddo
+
+case ('input.jgr')
+ call define_bc_jgr (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo)
+
+case ('input.sphere')
+ call define_bc_sphere (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo)
+
+case ('input.3Dpunch')
+ call define_bc_3Dpunch (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo)
+
+case ('input.2Dpunch')
+ call define_bc_2Dpunch (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo)
+
+case ('input.2Dpunch_vert')
+ call define_bc_2Dpunch_vert (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo)
+
+case ('input.folding')
+ call define_bc_folding (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo)
+
+case ('input.pipo')
+ call define_bc_pipo (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo)
+
+case ('input.ritske')
+ call define_bc_ritske (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo)
+
+case ('input.ritske_isurf')
+ call define_bc_ritske_isurf (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo)
+
+case ('input.subduction')
+ call define_bc_subduction (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo)
+
+case ('input.riedel')
+ call define_bc_riedel (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo)
+
+case ('input.parallipipede')
+ call define_bc_parallipipede (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo)
+
+case default
+ if (iproc.eq.0) print *,params%infile
+ call stop_run (' pb with input file in define_bc.f90')
+
+end select
+
+! here we need to impose that the bc satisfy the bad faces constrains
+
+do ie=1,osolve%nface
+ do j=1,4
+ jp=1+mod(j,4)
+ i=osolve%iface(j+4,ie)
+ if (osolve%kfix((i-1)*3+1).eq.1) osolve%u(i)=(osolve%u(osolve%iface(j,ie))+osolve%u(osolve%iface(jp,ie)))/2.d0
+ if (osolve%kfix((i-1)*3+2).eq.1) osolve%v(i)=(osolve%v(osolve%iface(j,ie))+osolve%v(osolve%iface(jp,ie)))/2.d0
+ if (osolve%kfix((i-1)*3+3).eq.1) osolve%w(i)=(osolve%w(osolve%iface(j,ie))+osolve%w(osolve%iface(jp,ie)))/2.d0
+ if (osolve%kfixt(i).eq.1) osolve%temp(i)=(osolve%temp(osolve%iface(j,ie))+osolve%temp(osolve%iface(jp,ie)))/2.d0
+ enddo
+i=osolve%iface(9,ie)
+if (osolve%kfix((i-1)*3+1).eq.1) osolve%u(i)=sum(osolve%u(osolve%iface(1:4,ie)))/4.d0
+if (osolve%kfix((i-1)*3+2).eq.1) osolve%v(i)=sum(osolve%v(osolve%iface(1:4,ie)))/4.d0
+if (osolve%kfix((i-1)*3+3).eq.1) osolve%w(i)=sum(osolve%w(osolve%iface(1:4,ie)))/4.d0
+if (osolve%kfixt(i).eq.1) osolve%temp(i)=sum(osolve%temp(osolve%iface(1:4,ie)))/4.d0
+enddo
+
+if (params%debug==2) call output_bc (osolve)
+
+return
+end
+!------------------------------------------------------------------------------|
+!------------------------------------------------------------------------------|
diff --git a/src/input.2011-01-21.txt b/src/input.2011-01-21.txt
new file mode 100644
index 0000000000000000000000000000000000000000..380aa4a4f080a2bf5121470b311b73a5ef7e0e3c
--- /dev/null
+++ b/src/input.2011-01-21.txt
@@ -0,0 +1,773 @@
+!------------------------------------------------------------------------------|
+!------------------------------------------------------------------------------|
+! |
+! ||===\\ |
+! || \\ |
+! || || //==\\ || || //==|| ||/==\\ |
+! || || || || || || || || || || |
+! || // || || || || || || || |
+! ||===// \\==// \\==\\ \\==\\ || |
+! |
+! 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_0000.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.07d-1
+
+[int]nstep is the number of time steps
+
+ nstep = 50
+
+[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 = -100
+
+[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 = 30
+
+[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 = 6
+
+[int]levelmax_oct is maximum level of octree discretization
+
+ levelmax_oct = 7
+
+[double]vex is the vertical exaggeration scaling factor that allows for
+variable aspect ratio elements. The value should range between 0-1, which
+will be multiplied by the vertical scale of the model domain. For example, a
+value of 0.25 would correspond to a model of dimensions 1x1x0.25, where the
+elements are 1/4 as tall as they are wide.
+
+ vex = 0.2d0
+
+[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=-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 = 7
+
+[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=1
+ npmax=4
+
+
+FEM + DIVFEM + MUMPS
+_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
+
+[int]matrule determines the rule used to assign material properties to
+elements that are cut by a surface. The default value (0) will use divFEM to
+determine the element properties. A value of 1 will assign the properties of
+the material that occupies the majority of the element, whereas a value of 2
+will assign the properties of the volumetric minority of the element. NOTE:
+divFEM is always used for any element containing void material.
+
+ matrule=1
+
+[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 = .026564d0
+
+ tempscale=100.d0
+
+
+ISOSTASY
+_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
+
+ isostasy = F
+
+ flexure = F
+
+ isobc = F
+
+MATERIALS
+_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
+
+[int]nmat is number of materials
+
+ nmat=2
+
+[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
+
+[bool]bulkvisc determines whether the user will input a material penalty value
+that scales with the shear viscosity or whether an independent bulk viscosity
+will be listed. If true, the values listed for [dp]penaltyi below should be
+bulk viscosity values for each material. Otherwise, the listed penalty value
+will be multiplied by the shear viscosity to enforce the incompressibility
+condition.
+
+ bulkvisc = F
+
+[bool]init_e2d determines whether an initial e2d value will be used to put
+plastic materials on yield in the first time step, or whether the standard
+initial material viscosity will be used. If true, an initial e2d value will be
+used to ensure plastic materials behave as desired.
+
+ init_e2d = T
+
+[dp]e2d0 is the initial strain rate value to use if init_e2d is true. This
+value will be used in only the first nonlinear iteration of the first time
+step and should be in scaled DOUAR units.
+
+ e2d0 = 5.d-4
+
+[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 strength
+- DPI-IV-> Drucker-Prager type of yield criterion
+ DPI -> the yield locus passes through the inner apices of the Mohr-Coulomb
+ hexagon
+ DPII -> the yield locus passes through the outer apices of the Mohr-Coulomb
+ hexagon
+ DPIII -> The average case, intermediate to DPI and DPII
+ DPIV -> Plane strain Drucker-Prager formulation
+ -> [dp]plasticity_1st_param is the angle phi
+ -> [dp]plasticity_2nd_param is the cohesion c
+ -> [dp]plasticity_3rd_param is the accumulated strain minimum, beyond
+ which phi is modified
+ -> [dp]plasticity_4th_param is the accumulated strain maximum, beyond
+ which phi is set to the value in plasticity_5th_param
+ -> [dp]plasticity_5th_param is the final phi value, if modified
+- DPV -> Drucker-Prager type of yield criterion that does not calculate the
+ values of alpha and k from phi and c. They are specified here. Note
+ that strain softening cannot be used for this DP type.
+ -> [dp]plasticity_1st_param is alpha
+ -> [dp]plasticity_2nd_param is k
+- DPVI -> Drucker-Prager type of yield criterion that uses the equation below
+ sigma_y = pressure * tan(phi) + c
+ -> [dp]plasticity_1st_param is the angle phi
+ -> [dp]plasticity_2nd_param is the cohesion c
+ -> [dp]plasticity_3rd_param is the accumulated strain minimum, beyond
+ which phi is modified
+ -> [dp]plasticity_4th_param is the accumulated strain maximum, beyond
+ which phi is set to the value in plasticity_5th_param
+ -> [dp]plasticity_5th_param is the final phi value, if modified
+- DPVII -> Drucker-Prager type of yield criterion that uses the same 2D
+ formulation that is used in Sopale
+ (sigma_y = pressure*sin(phi) + c*cos(phi))
+ -> [dp]plasticity_1st_param is the angle phi
+ -> [dp]plasticity_2nd_param is the cohesion c
+ -> [dp]plasticity_3rd_param is the accumulated strain minimum, beyond
+ which phi is modified
+ -> [dp]plasticity_4th_param is the accumulated strain maximum, beyond
+ which phi is set to the value in plasticity_5th_param
+ -> [dp]plasticity_5th_param is the final phi value, if modified
+- MC -> Mohr-Coulomb type of yield criterion
+ -> [dp]plasticity_1st_param is the angle phi
+ -> [dp]plasticity_2nd_param is the cohesion c
+ -> [dp]plasticity_3rd_param is the accumulated strain minimum, beyond
+ which phi is modified
+ -> [dp]plasticity_4th_param is the accumulated strain maximum, beyond
+ which phi is set to the value in plasticity_5th_param
+ -> [dp]plasticity_5th_param is the final phi value, if modified
+
+ 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 = -1.d0
+ viscosity1 = 1.d5
+ penalty1 = 1.d8
+ expon1 = 1.d0
+ diffusivity1 = 5.98d-4
+ heat1 = 0.d0
+ activationenergy1 = 0.d0
+ plasticity_type1 = DPVII
+ plasticity_1st_param1 = 20.d0
+ plasticity_2nd_param1 = 0.d0
+ plasticity_3rd_param1 = 0.5d0
+ plasticity_4th_param1 = 1.5d0
+ plasticity_5th_param1 = 20.d0
+
+ density2 = -.935d0
+ viscosity2 = 1.d-3
+ penalty2 = 1.d8
+ expon2 = 1.d0
+ diffusivity2 = 1.20d-3
+ heat2 = 0.d0
+ activationenergy2 = 0.d0
+ plasticity_type2 = No
+
+[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.d7
+
+
+SURFACES
+_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-
+
+[int]ns is number of surfaces to track
+
+ ns=2
+
+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
+ 12 for a rectangular emboss with defined margin slope
+ -> 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
+ -> surface_param_07 is the slope in degrees
+ - [int]leveloct is the level at which the octree will be refined in the vicinity
+ of the surface.
+
+ levelt1 = 8
+ itype1 = 0
+ surface_type_1 = 12
+ rand1 = F
+ surface_param_01_1 = 0.01289d0
+ surface_param_02_1 = -0.5d0
+ surface_param_03_1 = 0.5d0
+ surface_param_04_1 = -1.d0
+ surface_param_05_1 = 2.d0
+ surface_param_06_1 = -0.01289d0
+ surface_param_07_1 = 1.47676d0
+ material1 = 1
+ activation_time_1 = -1.
+ leveloct1 = 7
+ stretch1 = 1.5d0
+ anglemax1 = 180.d0
+ criterion1 = 1
+ anglemaxoctree1 = 180.d0
+ spread_surface_points1 = 1
+
+ levelt2 = 8
+ itype2 = 0
+ surface_type_2 = 12
+ rand2 = F
+ surface_param_01_2 = 0.d0
+ surface_param_02_2 = 0.0625d0
+ surface_param_03_2 = 0.9375d0
+ surface_param_04_2 = -1.d0
+ surface_param_05_2 = 2.d0
+ surface_param_06_2 = -0.01289d0
+ surface_param_07_2 = 87.5d0
+ material2 = 2
+ activation_time_2 = -1.
+ leveloct2 = 7
+ stretch2 = 1.5d0
+ anglemax2 = 180.d0
+ criterion2 = 1
+ anglemaxoctree2 = 180.d0
+ spread_surface_points2 = 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 = 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.d0
+ box1z1=.03d0
+ box1level=7
+
+ box2x0=0.d0
+ box2x1=.75d0
+ box2y0=.4d0
+ box2y1=.6d0
+ box2z0=.05d0
+ box2z1=.25d0
+ box2level=7
+
+ 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=.0428d0
+
+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=-600.d0
+ velocity_scale=50.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