Sopale documentation

USER GUIDE for SOPALE1 (Serial OPALE)

Version: code-mar17-04

INPUT FILES

   SOPALE1_i
   sopale1_plus_i
   surfaceproc_i
   matlab_i (optional)
   
 The following files haven't been documented yet (most are for Philippe's
 debugging) but SOPALE still needs them, so just copy them from an existing
 run:

   SOPALE1_interactive_graphics_i
   analyze_wedge_growth_i
   check_picard_i
   check_picard_show



File:SOPALE1_i 

   read(8,*)template_i
        template_i = number of lines to skip
                     used for reading header= model templates listed
                     at the beginning of the file
  
   do i=1,template_i
      read(8,'(a1)')modelkey
   enddo
       This part skips over the lines which describe the model templates

  read(8,*)nmodels
        nmodels = number of models

  read(8,*)(modelsn(i),i=1,nmodels)
        list of model numbers to run
      
  For each model, we have the following input:
      
      read(8,'(a1)')modelkey
           modelkey = one-character identifier for the model        
      
      read(8,*)amodel, benchmark
           model number for this model, benchmark number
           
      read(8,*)wres,rres, align_time
           wres = number of timesteps between writing of 2 restart files
           rres = determines whether to read restart file
                  rres = 0 - don't read restart file
                  rres = 1 - read restart file
           align_time (only has effect when rres = 1)
                      = 0  on restart, set start time to zero and timestep to 1
		           This is appropriate when you're using a restart file
			   which just establishes the state you want.
                      = 1  on restart, set start time to time in restart file
                           and set timestep to timestep in restart file.
			   This is good for restarting an interrupted run.
           
      read(8,*)nsave1
           number of saves for Eulerian grid 
           (Maximum number of saves is 200)
      
      read(8,*)(isave_1(i),i=1,nsave1)
           [only needed if nsave1 > 0 ]
           timesteps at which to save the Eulerian grid
               
      read(8,*)nsave2
           number of saves for Lagrangian grid
           (Maximum number of saves is 200)
           
      read(8,*)(isave_2(i),i=1,nsave2)
           [only needed if nsave2 > 0 ]
           timesteps at which to save the Lagrangian grid
      
      read(8,'(a14)')nameout1
            14-character prefix for the name of the output files
            NOTE: there must be NO MORE and NO LESS than 14 characters,
            and the last two characters must be "/g".
            
            Bonny's note:  Philippe was assuming that the prefix would
            always be './Sopale/rXX/g' where XX is 01, 02, 03, etc.
            Because of the way Philippe wrote the code, you cannot use a
            prefix of more than or less than 14 characters. Using blanks
            before or after a shorter prefix will give problems with 
            embedded blanks in filenames.  This means that if you
            want to keep the input files and the output file in the
            same directory, you would have to use the prefix
            '././././././/g'.

      read(8,*)time_nstep,dt,outunt
            time_nstep = number of timesteps in this model run
            dt         = timestep length
            outunt     = Fortran unit number for output, e.g. 60
            
      read(8,*)bigu,bigt
            bigu = ~1d7 multipliers for penalizing boundary velocity
            bigt = ~1d7 multiplier for penalizing the temperature boundary
                   conditions
      
      read(8,*)length,height,nx1,ny1,betaxg,betayg,length2,height2,nx2,ny2
           length = length of Eulerian rectangle
           height = height of Eulerian rectangle
           nx1 = number of rows of Eulerian grid
           ny1 = number of columns of Eulerian grid
           betaxg = 1 for regular grid (heterogeneous grids are other values)
           betayg = 1 for regular grid (heterogeneous grids are other values)
           length2 = length of Lagrangian rectangle
           height2 = height of Lagrangian rectangle
           nx2 = number of rows of Lagrangian grid
           ny2 = number of columns of Lagrangian grid
 
           NOTE: surface is at "height", not 0
                 nodes are numbered starting at top surface
                 x=0 is at left edge 
            
             
     
      read(8,*)npert
           npert = number of interfaces to perturb for instability analysis
                 = 0 for sandbox analysis
      
      read(8,*)(ipert(i),i=1,npert)
           [only required if npert > 0]
      
      read(8,*)(lpert(i),i=1,npert)
           [only required if npert > 0]
      
      read(8,*)(lpert0(i),i=1,npert)
           [only required if npert > 0]
      
      read(8,*)(apert(i),i=1,npert)
           [only required if npert > 0]
      
      read(8,*)nfpert
           nfpert = perturbation of forces
                  = 0 for our models
                  
      read(8,*)afpert,flpert0,flpert,fipert,xyfpert
           [only required if nfpert > 0]
            
      read(8,*)iformul

      read(8,*)(hists(i),i=1,11)

      read(8,*)ang_gravity,gravity,pcohes,visminmax(1),visminmax(2),psmooth,
           p_bdry_layer, splasmin, vismin_bdry_layer, edotmax

           ang_gravity = angle of gravity in degrees
			 The convention is that ang_gravity is positive
			 when measured anti-clockwise from the positive
			 X direction. Gravity directed in positive
			 X direction has ang_gravity=0.
           gravity = value of gravity
                     assumed to act toward decreasing Z
           pcohes  = lower bound for pressure (pressures below this 
	             are set to pcohes)
           visminmax(1) = minimum viscosity limit (any viscosities
                          below this are set to visminmax(1))
           visminmax(2) = maximum viscosity limit (any viscosities
                          above this are set to visminmax(2))
           psmooth = 1 for pressure smoothing
                   = anything else, no pressure smoothing

           p_bdry_layer = minimum pressure for setting minimum splas and
                          viscosity; i.e. if the pressure is less than 
                          p_bdry_layer, then the splas and viscosity are
                          kept above the minimums specified in splasmin
                          and vismin_bdry_layer.
           
           splasmin  = minimum splas limit for areas where epress(i) is less
                       than p_bdry_layer
           
           vismin_bdry_layer  = minimum viscosity limit for areas where 
                                epress(i) is less than p_bdry_layer.

           edotmax = maximum value of the second invariant of the
                     strain rate

                          **** NOTE on using edotmax ****
                     During iterations when the frictional sediment 
                     near the surface becomes unstable and starts to slump,
                     the model viscosity in this region is progressively
                     reduced during iterations in an attempt to put the
                     material on plastic yield.  This numerical softening
                     leads to run-away in the viscosities in the 
                     subroutine viscosities1, where viscosity is set to
                     sleast/edot, where sleast is either the yield stress
                     for a plastic material or the viscous stress. Increasing
                     velocities (e.g. slumping) lead to increases in edot
                     and then to decreases in viscosity which then results
                     in a further increase in velocity and the cycle spirals
                     down out of control with decreasing viscosities.
                     This does not occur for viscous materials because they
                     always are assigned their true material viscosity, not
                     an artificial one that puts them on plastic yield.
                       The use of edotmax limits the maximum strain rate
                     and controls the velocity of slumping.  It does not
                     guarantee good convergence during iterations but it
                     allows controlled slumping so that the model can be
                     continued through a slumping event.  All other parts of
                     the model will be nicely converged.
                       The way to run this is to choose an edotmax that is
                     larger than any of the normal strain rates in the
                     model, so it won't affect them. Set the max number of
                     iterations per timestep to a relatively low value
                     so you don't have a large number of timesteps with many
                     iterations, and set force_continue to be 'on' so 
                     that the model will continue to timestep even though
                     the convergence is not as good as the criterion you have
                     set during slumping events.
      

      read(8,*)color_number
           color_number = number of materials whose mechanical properties 
                          are to be defined
      
      read(8,*)mixhyp,refsrate,redphi,nredphi, plas_fluid_method

             mixhyp = mixing hypothesis for viscous material properties
		      (viscous rheologies)
                    = 1  viscous rheologies act in serial
                         NOT CODED
		    = 2  viscous rheologies act in parallel.
			 Each viscous rheology is assumed to have the
			 same second invariant of the strain rate (they
			 act in parallel). Stresses are determined
			 from their respective effective viscosities X
			 strain rate.  Total second invariant of stress
			 is sum of individual stresses.
			 The total viscous stress is compared with the
			 plastic yield stress to determine whether
			 flow is viscous or plastic.
		    = 3  minimum dissipation
			 Viscous rheologies combine in a way that
			 minimizes dissipation
			 NOT CODED

                    Use mixhyp = 2

             refsrate = reference strain rate to start viscosity computation
                        for linear creep.

			NOTE: for certain choices of the parameters in each
			of the nonlinear creep models (gscreep,nlcreep,pei)
			the calculation for the initial viscosities may be
			insensitive to the value of refsrate. Sigma1, sigma2,
			and sigma3 contain a multiplier edot and viscos is
			calculated by dividing the total stress by 2*edot, so
			that viscos is insensitive to the choice of edot under
			these circumstances.
			See the code in subroutine viscosities1, in the branch 
			where notinitialize.ne.1.

             NOTE: special method for setting mechanical material properties,
             used primarily to create a bottom boundary layer. This sets
             the material properties for the mechanical calculations only,
             it does NOT change the actual color (material) number of the
             Eulerian elements.  Set nredphi > total number of Eulerian
             elements to suppress this method.
             
             redphi = mechanical material number for all elements
                      greater than nredphi
             nredphi = element number.  Eulerian elements greater
                       that this number will use the properties of
                       material redphi for mechanical calculations.
       		       Note: elements are numbered by row, so this
		       can be used to set a basal boundary layer. The
		       mechanical material properties overrides
		       those specified through the normal input
		       process (Quick and dirty method for setting
		       all elements to material redphi - use redphi=
		       0)
             
	     CB Comments: redphi and nredphi are used to define
	     special conditions for the mechanical Eulerian element
	     materials. If the element number is greater than
	     nredphi, it will be assigned the mechanical material
	     properties of material number redphi for mechanical 
             calculations.  Note that the actual color number (material
             number) of the Eulerian element is not overwritten, so 
             you will not see the boundary layer material redphi in plots
             of material colors.

	     Use of redphi and nredphi overrides other methods used
	     to assign mechanical properties. [Example of use is to
	     define a basal boundary layer in sandbox type
	     experiments such that the boundary layer is redefined
	     every timestep according to the Eulerian grid and
	     therefore in a manner that overrides material advection
	     and regridding from the Lagrangian grid.]
	     
	     nredphi should be used with care, remember that it
	     assigns one mechanical material property (that of
	     redphi) to the Eulerian elements. The effect of strain
	     dependent material properties should also be considered.
	     To avoid resetting element properties in this way, set
	     nredphi.gt.max number of Eulerian mechanical elements.

	     plas_fluid_method  = method for calculating plastic yielding
	                          for a frictional plastic
	                          incompressible material using a
	                          Drucker-Prager yield, for
	                          plane-strain and including various
	                          options for the effects of pore
	                          fluid pressure (see notes below).
	     
	     plas_fluid_method not specified in input
	        or
	     plas_fluid_method = 0
	         Defaults to old behaviour in which phi values (see
	         below) are interpreted as effective phi's modified
	         by pore fluid pressure.
	         
	         In subroutine viscosities1, yielding is given by
	            splas = epressure * dsin(phi) + coh

	     plas_fluid_method = 1
	         Assume lambda_pore_prs = fluid pressure/solid pressure
	                                                 (mean stress)
	           splas = epressure*(1-lambda_pore_prs)*dsin(phi) + coh
                 (small error in last term)

	     plas_fluid_method = 2
	         Assume lambda_pore_prs = fluid pressure/solid pressure
	                                                 (mean stress)
	           splas = epressure*(1-lambda_pore_prs)*dsin(phi)
	                   + coh*dcos(phi)
                 (small error corrected)

	     plas_fluid_method = 3  (see warning below)
                 Assume lambda_pore_prs = fluid pressure/overburden weight
                 (in code, overburden weight is called plitho)
                 
	           splas = (epressure - lambda_pore_prs*plitho)*dsin(phi) + coh
                 (small error in last term)

	     plas_fluid_method = 4  (see warning below)
                 Assume lambda_pore_prs = fluid pressure/overburden weight
                 (in code, overburden weight is called plitho)
                 
	           splas = (epressure - lambda_pore_prs*plitho)*dsin(phi)
	                   + coh*dcos(phi)
                 (small error corrected)

	     plas_fluid_method = 5 
                 The effect of hydrostatic fluid pressures
                 on epressure is taken into account.  Method assumes
                 that the fluid pressure is equal to the hydrostatic
                 fluid pressure everywhere.  Because fluid pressure is
                 hydrostatic, lambda_pore_prs defaults to 1.
		 (in code, hydrostatic weight is called phydro)

	          splas = (epressure - lambda_pore_prs*phydro)*dsin(phi)
	                   + coh*dcos(phi)
	         
	         When you use plas_fluid_method=5, lambda_pore_prs1
	         and lambda_pore_prs2 must be set to 1 in SOPALE1_i.
	         	         
	         NOTE: The hydrostatic pressure (phydro) is applied to
	         ALL materials in the model.  User beware!


             NOTE: Methods 1 and 3 have errors in the cohesion term.
             These methods are available for comparison with models
             that used the old behaviour.  Methods 2 and 4 have the
             correct cohesion term.
             
             WARNING!  As of August 2002, models using 
             plas_fluid_method = 3 and plas_fluid_method = 4 will not
             converge because epressure can be less than
             lambda_pore_prs*plitho, resulting in a negative splas. 
             A fix is still under consideration.
             [mod BL 2002/09/18]
	     ** UPDATE **   Using lambda_blend1 and lambda_blend2 to
	     gradually introduce the lambda_pore_prs seems to improve the
	     chances of the model getting past the first timestep.  Also,
	     values of lambda_pore_prs greater than 0.4 seem to cause
	     problems. 
             [/mod BL 2002/09/18]
             

       For each material (i=1,color_number), we have the following input:
                                     
         read(8,*)density(i)
              density(i) = density of this material
      
	 read(8,*)bulkv(i)
	      bulkv(i) = bulk viscosity

         [mod 2002/07/23]
	     NOTE:
             Plasticity yield criterion updated July 23, 2002 to include
             fluid pore pressure effects explicitly.
             Plastic yielding uses the Drucker-Prager criterion
             
                       1/2
                  (J ')    = p(1-lambda) X sin(phi) + c
                    2

             where J ' is the second invariant of the deviatoric stress
                    2   
                   p   is the pressure in the solid
                   phi is the internal angle of friction
                   c   is the cohesion
                   lambda is the pore pressure ratio (p      / p     )
                                                       fluid    solid
                                                       
	     In subroutine viscosities1, this appears as
	     
	       splas = epressure*(1.0-lambda_pore_prs)*dsin(phi) + coh
	     
	     (Old version was  splas=epressure*dsin(phi) + coh)
	     
	     phi, lambda_pore_prs, and coh are material properties
	     specified in input for each material (see below,
	     description of next two lines of input).  phi,
	     lambda_pore_prs and coh can be specified to vary with
	     the second invariant of the deviatoric strain, here
	     termed "strain".
	       
	 read(8,*)phi1(i),phi2(i),phie1(i),phie2(i),
	          lambda_pore_prs1(i), lambda_pore_prs2(i),
	          lambda_pore_prse1(i), lambda_pore_prse2(i)
	          
	      phi1(i) = first angle of friction in degrees
	      phi2(i) = second angle of friction in degrees
	      phie1(i) = strain level 1
	      phie2(i) = strain level 2
	      NOTE: must have phie1(i) < phie2(i)
	      
	      Diagram:
	      
	                              ________________ phi2
	                             /|
	                            / |
	                           /  |
	                          /   |
	      phi1 ______________/    phie2 (strain level 2)
	                         |
	                        phie1
	                        (strain level1)    
	      
	      lambda_pore_prs1(i) = first pore pressure ratio
	      lambda_pore_prs2(i) = second pore pressure ratio
	      lambda_pore_prse1(i) = strain level 1
	      lambda_pore_prse2(i) = strain level 2
	      NOTE: must have lambda_pore_prse1(i) < lambda_pore_prse2(i)
	      
	      
	                              ________________ lambda_pore_prs2
	                             /|
	                            / |
	                           /  |
	                          /   |
      lambda_pore_prs1 __________/    lambda_pore_prse2
	                         |
	                        lambda_pore_prse1     

         
         Use of phi and lambda_pore_prs
         
         Previous use varied phi to take account of lambda_pore_prs
         implicitly, i.e. by choosing phi_effective such that
            sin(phi_effective) = (1.0-lambda_pore_prs)*sin(phi)
         This method can still be used in which case all inputs for
         labmda_pore_prs should be set to 0 and the input values of
         phi are interpreted as effective internal angles of
         friction.
         
         Alternatively, the input phi's are true, material
         properties, e.g. compatible with Byerlee's Law (30-40
         degrees) and the labmda_pore_prs values are effective pore
         pressure ratios associated with each material.  Note fluid
         pressures are not calculated separately/dynamically in this
         verison.  Lambda_pore_prs is specified as a material property
         and should be chosen to reflect the fluid pressure
         environment envisioned for this material.

         lambda_pore_prs =
             0.0      dry rock environment, no fluid pressure
             0.3-0.4  hydrostatic fluid pressure - depends on the
                      density distribution of the model material
             close to 1.0 (but not equal to 1.0)
                      lithostatic conditions, strongly overpressured
                      conditions
         
	 Note that the theory used here is the most simple form.  It
	 hids the underlying problem of the best measure of the solid
	 pressure in lambda_pore_prs.  The solid pressure may not be
	 lithostatic.  This is likely the case in frictional-plastic
	 materials (see Petrini and Podladchikov, J.Met.Geol.v.18,
	 p.67-77,2000).

      [/mod 2002/07/23]
      
      [mod 2002/08/01]

 NOTES CONCERNING
 CHOICE OF FRICTIONAL PLASTIC FAILURE CRITERIA IN SOPALE-
  METHODS 1 AND 2 FOR INCLUDING PORE FLUID PRESSURE EFFECTS USING
  THE HUBBERT-RUBEY PORE FLUID PRESSURE RATIO
 
  Derivation is 2D and then generalised to Plane Strain  

  Frictional Plastic Failure Criteria Including Effects of Pore Fluid
   Pressure

   Note in this section the names of variables/parameters are
   not the same as those used elsewhere in the guide. The 
   approx conversion is
     phi = f
     lambda_pore_prs = L 
     coh = C   
     stress = S (various components)

1)  Definition of Hubbert-Rubey Pore Fluid Pressure Ratio

  As Glen points out Hubbert and Rubey have defined L (or lambda)
   as 
     L1 = Pf / Sn
 
   and as
     L2 = Pf / Sz

   where Pf is the pore fluid pressure, Sn is the normal stress
   component acrossa potential failure plane and Sz is the 
   weight of the overburden ( equivalent to the lithostatic weight
   rho x g x z in a uniform medium).

  Here we use the second definition in the same manner as Davis et al
  1983. They define the vertical normal traction to be solely the
  effect of the lithostatic overburden and include a hydrostatic
  component for submarine problems.

   To avoid confusion it is probably best to refer to Sz as the 
  weight of the overburden. The confusing term is the 'lithostatic
  pressure'...which is sometimes used to mean the weight of the
  overburden and sometimes used to mean the solid pressure.    
    
   The stress convention is that compressive stresses are positive
  ( this is not the normal convention...but it is also used by 
   Davis et al and Jaeger and Cook.)
 
2) Definition of Effective Stress
     Se(1) = S(1) - Pf,    Se(3) = S(3) - Pf ,

   where Se(1), S(1) are the maximum effective principal stress and
   the maximum principal stress and the Se(3) and S(3) are the 
   corresponding minimum principal stresses.

3) The Coulomb Criterion (Jaeger and Cook, page 95)

   Shear failure in rock occurs when

    Mod( T ) = C + S x Tan (f)                                (1
   where S and T are the normal and shear stresses across the failure
   plane, C is a constant, the cohesion and f is the internal
   angle of friction ( normally called phi).

   This can also be expressed in terms of the Mohr circle
   representation.

    1/2 x (S(1)-S(3)) = 1/2 x(S(1)+S(3))x Sin(f) + C x Cos(f) (2

4) Effect of Fluid Pressure on the Coulomb Failure Criterion 
 
   The Mohr envelope for the Coulomb criterion is modified by
    pore fluid pressure to

      Mod (T) = (S - Pf) x Tan (f) + C
   or 
      Mod (T) = Se x Tan (f) + C   
   where Se is the effective normal stress on the failure plane.

   Similarly the Mohr circle representation can be written in
   terms of effective stresses

   1/2 x( Se(1)-Se(3)) = 1/2 x (Se(1)+Se(3))x Sin(f) + C x Cos(f) (3  

 5) Mean Solid Pressure

    the mean pressure in the solid is

       P = 1/2 x (S(1)+S(3))

    Therefore  P = 1/2 x (Se(1)+Se(3)) + Pf                  (4
    Substituting for Se(1)+Se(3) from 4 into 3

    1/2 x(Se(1)-Se(3)) = (P - Pf) x Sin(f) + C x Cos(f)       (5
    Because 
    Se(1)-Se(3) = S(1)-S(3)
    5 can be written

    1/2 x(S(1)-S(3)) = (P - L x Sz) x Sin(f) + C x Cos(f)

    This is the two dimensional form of the Coulomb failure 
     criterion taking pore fluid pressure into account through
     the Hubbert-Rubey pore fluid pressure ratio. 

 6) Failure Criterion in Plane Strain in Terms of the Stress 
    Invariants

     If it is assumed that in plane strain that the intermediate
    principal stress acts in the direction normal to the plane
    containing the other principal stresses and has a value
    approximately equal to the solid pressure then both the 
    second invariant of the stress and the second invariant
    of the deviatoric stress both reduce to

     Sqrt (I2) = Sqrt (J2D) = 1/2 x (S(1)-S(3))

    where Sqrt is the square root.
     ASIDE
    This is because for S(2) = P it follows that 
     S(1)-S(3) = 2 x(S(2)-S(3)) = 2 x (S(1)-S(2)) and 
      therefore 
      Sqrt (I2)= 1/6 x [(S(1)-S(3))^2 + (S(3)-S(1))^2 
                             + (S(1)-S(2))^2]
      becomes  
               = 1/6x[ (1/4  + 1 + 1/4)(S(1)-S(3))]
               = 1/2 x(S(1)-S(3))....as above
     Note this is only correct if the intermediate principal
      stress is out of the plane and equal to the pressure.

     It neeeds to be shown that this is also true if yielding
     is written in terms of the second invariant of the deviatoric
     stress.

7) 'Mohr-Coulomb' Failure Law
     Montesi and Zuber give the M-C failure law in terms
      of the invariants as
 
       I2 = I1 x Sin(f) + C x Cos(f) their eqn 18
     This is close to the form given above but differs in
     sqrt(I2) vs I2  and I1 vs P. The eqn is not dimensionally
     correct with I2 .They may choose definitions of I2 as sqrt(I2)
     and I1 = I1/3 = P 

     END ASIDE   

    Under these circumstances the failure criterion can be 
    written as 

      Sqrt (J2) = Sqrt (J2D) = (P - L x Sz) x Sin(f) + C x Cos(f)

    where P = 1/3 (I1) = 1/3 x ((S(1) + S(2) + S(3)) 

     The above is the failure criterion when Plas_Fail_ Method = 2
    and appears in subroutine viscosities as 

    splas = (epressure - lambda_pore_prs x plitho) x dsin(phi)
        + coh 

     Note the small error that coh  should be coh x dcos(phi)

    END NOTES
      [/mod 2002/08/01]
      	      
	 read(8,*)coh1(i),coh2(i),cohe1(i),cohe2(i)
	      coh1(i) = first cohesion
	      coh2(i) = second cohesion
	      cohe1(i) = strain level 1
	      cohe2(i) = strain level 2
	      NOTE: must have cohe1(i) < cohe2(i)
	      
	      Diagram:
	      
	                              ________________ coh2
	                             /|
	                            / |
	                           /  |
	                          /   |
	      coh1 ______________/    cohe2
	                         |
	                        cohe1     
	      
	      	      
	 [GRAIN SIZE SENSITIVE CREEP]
	 read(8,*)gscreep_mu(i),gscreep_m(i),gscreep_T(i),
	        gscreep_strain1(i), gscreep_strain2(i),gscreep_a0(i),
	        gscreep_a1(i), gscreep_a2(i), gscreep_updown_T
	        
	  [1] gscreep_mu(i) = viscosity at level of grainsize a0

                             NOTE: See Philippe's note on viscosities
			     in SOPALE (there is also a link to this from
			     the sopale documentation web page).

	  [2] gscreep_m(i) = power of grainsize dependency
	                     [m* in formula below]
	  [3] gscreep_T(i) = temperature at which mu is (1/e) x that viscosity
	                     at a0 and T
	                     [T* in formula below]
	  [4] gscreep_strain1(i) = S1* (see fig below)
	  [5] gscreep_strain2(i) = S2* (see figure below)
	      NOTE: Always use S1* < S2*
	      
	  [6] gscreep_a0(i) = parameter a0* [ a*  in formula below] 
	                                       0
	  [7] gscreep_a1(i) =  a1* (see fig below)
	  [8] gscreep_a2(i) =  a2* (see fig below)
	  [9] gscreep_updown_T 
	   
	                               
	                               [2]
	              d                m*         d
	          mu(T  ) = mu* (a/a* )   exp (- T  /T*)
	                            0                [3]  
	                           [6]                
                 d = dynamical



	    [8] a1* _ _ _ _ _ _ _ _ _ ______________
	                             /|
	                            / |
	                           /  |
	                          /   |
	    [7] a2* _____________/    |
	                         |    |
	                         |    | 	          
	                        S1*   S2*
	                        [4]   [5]
	          
	                       d
	     For case  exp (- T  /T*)   gscreep_updown_T = +1
	     
	                          d
	     For case  exp (- T*/T  )   gscreep_updown_T = -1
	          
	  
	  
	  To suppress gs_creep in parallel, take gs_creep_mu = 0
	     
	 ******** insert note about Rigid-Plastic calculation here?? ****
	  
	     
	     
	              
	 read(8,*)nlcreep_eps(i),nlcreep_sc1(i),nlcreep_sc2(i),
	          nlcreep_n(i),nlcreep_T(i),nlcreep_strain1(i),
	          nlcreep_strain2(i),nlcreep_updown_T
	          
	    [1] nlcreep_eps(i)
	    [2] nlcreep_sc1(i)
	    [3] nlcreep_sc2(i)
            [4] nlcreep_n(i)
            [5] nlcreep_T(i)
            [6] nlcreep_strain1(i)
            [7] nlcreep_strain2(i)
            [8] nlcreep_updown_T
	        
	      See formula and diagram below...
	        
	                             [4]      
	        .   .                n*             l
	      ( e / e* ) = (sigma/sc)    exp -(T*/ T  )
	             0             |          [5]
	            [1]           sc1=[2]          if updown= -1
	                          sc2=[3]         



        	[2]  sc1* _____________ 
                        	      |\
                        	      | \
                        	      |  \
                        	      |   \
        	[3]  sc2* _ _ _ _ _ _ | _ _\__________________
                        	      |    |
                        	      |    |
                        	      S1*  S2*
                        	      [6]  [7]


          PEIERL'S LAW POWER-LAW BREAKDOWN
	          
	 read(8,*)pei_eps(i),pei_Tn(i),pei_sc(i),pei_T2(i),
	          pei_sc(i),pei_T2(i),pei_strain1(i),pei_strain2(i) 
	 
	    [1] pei_eps(i) = reference strain rate
	    [2] pei_Tn(i)    [use 0 to suppress in parallel]
	    [3] pei_sc(i)
	    [4] pei_T2(i)
	    [5] pei_sc(i)
	    [6] pei_T2(i)
	    [7] pei_strain1(i)  - not used
	    [8] pei_strain2(i)  - not used
	 
	 
	  NOTE: E = epsilon
	  
	                           [2]
	         .   .             Tn*/T          [4]
	         E = E* (sigma/sc*)       exp ( - T2* / T ) 
	             [1]       [3]                        
	 
	 
	 read(8,*)texpand(i),tdensity(i),
                  melt_weaken(i), temp_solidus(i),
                  temp_melt_weakened(i), visc_melt_weakened(i)

	       texpand(i) = volume coefficient of thermal expansion
	       tdensity(i) = reference temperature for density 
			     (in absolute degrees)
	                     if positive, density is
	                        rho(T) = rho (1 - alpha(T - Tref))
	                                    0
                                [where Tref is tdensity(i)]
	                     if negative, density if
	                        rho(T) = - rho  alpha T
	                                      0
                melt_weaken(i) = 1 if material melt weakens
                               = 0 if no melt weakening
                temp_solidus(i) = temperature (degrees Kelvin) at which 
                                  melt weakening starts
                temp_melt_weakened(i) = temperature (degrees Kelvin) at
                                        which melt weakening is complete
                visc_melt_weakened(i) = linear viscosity of melt-weakened
                                        material

	 [END OF INPUT FOR EACH MATERIAL'S MECHANICAL PROPERTIES]
	 
  
      read(8,*)type_bc
  
      IF type_bc = 1 then

           type_bc = 1 for singularity model   --------o---------
                                                       S

      read(8,*)vleft,vright,xsleft,xsright
           vleft = velocity at far left of singularity
           vright = velocity at far right of singularity
           xsleft   }  [see diagram]  velocity changes linearly from vleft
           xsright  }                 to vright across the zone
           
               xsleft            xsright
                 |                 |
               --+--------o--------+---  
                 |        S        |
          

	This type of bc does not allow vertical velocities on the
	left or right boundaries; i.e. these are not roller boundaries.
	Use this type bc for models where the deformation does not
	reach to the ends of the Eulerian mesh.
	NOTE: this means there is no need to modify the Lagrangian
	grid that is outside the Eulerian mesh owing to vertical
	movements at the left and right boundaries.

	For cases where the left and/or right boundaries move vertically,
	use type_bc = 11.
      
              -------------------------------------------------
                THE TYPE_BC = 3 FAMILY OF BOUNDARY CONDITIONS
              -------------------------------------------------
             
        There are a number of boundary conditions based on type_bc=3,
        which have numbers starting with 3 (i.e. type_bc=30,
        type_bc=31, type_bc=312 and so on).  They are variations on
        a theme, and have some basic properties in common, as
        illustrated by this diagram:
            


	    +--------------------------------------+
            |-->                                <--|
            |--> vxleft_up           vxright_up <--|
            |-->                                <--|
            |                                      |
    yleft1 -+-                                    -+- yright1
            |                                      |
            |                                      |
    yleft2 -+-                                    -+- yright2
            |                                      |
         <--|                                      |-->
         <--| vxleft_down             vxright_down |-->
         <--|                                      |-->
         <--|                                      |-->
         <--|                                      |-->
         <--|             vybottom                 |-->
         <--|           ^ ^ ^ ^ ^ ^ ^ ^            |-->
            |        ^  | | | | | | | |  ^         |
            |     ^  |  | | | | | | | |  |  ^      |
            +----+---------------------------+-----+
                 |                           |
               xvyleft                     xvyright
               
        NOTE: The bottom boundary y velocity is set to vybottom
	      between xvyleft and xvyright, zero to the left of xvyleft,
	      and zero to the right of xvyright.  The variable
	      decaylvybottom sets the length for the bottom boundary y
	      velocity to decay to zero.
		In the Lagrangian grid below the Eulerian grid (i.e. where
	      the y coordinate of the grid is less than 0), the
	      y velocities are set to the bottom boundary y velocity.
	      For now (as of June 4, 2002) we are setting the x
	      velocities in this part of the Lagrangian grid to zero.
	      This means that the mechanical and thermal materials in 
	      in this area of the Lagrangian grid are advected
	      but not deformed.


	vxleft_up = x velocity at left boundary above yleft1
	vxleft_down = x velocity at left boundary below yleft2
              
	      In the area between yleft1 and yleft2, the
	      boundary condition makes a smooth transition between
	      vxleft_up and vxleft_down.
      
	yleft1 = vertical position at left boundary, above which
	         the x velocity is vxleft_up
	yleft2 = vertical position at left boundary, below which
	         the x velocity is yleft_down

	rollleft = 1  means roller (free slip) boundary condition on
		      left boundary
		 = 0  means no vertical component to boundary condition
		      on left boundary, i.e. pinned (no slip)

	vxright_up = x velocity at right boundary above yright1
	vxright_down = x velocity at right boundary below yright2
              
	      In the area between yright1 and yright2, the
	      boundary condition makes a smooth transition between
	      vxright_up and vxright_down.
      
	yright1 = vertical position at right boundary, above which
	         the x velocity is vxright_up
	yright2 = vertical position at right boundary, below which
	         the x velocity is yright_down

	rollright = 1  means roller (free slip) boundary condition on
                       right boundary
		  = 0  means no vertical component to boundary condition
		       on right boundary, i.e. pinned (no slip)
		      
        vytop0    = 1  to set a rigid top (i.e. no vertical
                       velocities on top boundary)
                  = 0  vertical velocities allowed on top boundary
        
        vxbottom0 = 1  for no slip on bottom (no horizontal
                       velocities on bottom)
                  = 0  horizontal velocites allowed on bottom

        vybottom = y velocity at bottom boundary   
        decaylvybottom = length scale in SI, m, over which 
                         vybottom (the bottom boundary y velocity)
                         decays to zero adjacent to xleft and
                         xright
        xvyleft = to the left of this point, bottom boundary y
                 velocity is zero
        xvyright = to the right of this point, bottom boundary y
                 velocity is zero

                  --------------------------------
		      
      IF type_bc =3 then
            
      read(8,*)vxleft_up,vxleft_down,yleft1,yleft2,rollleft,
        vxright_up,vxright_down,yright1,yright2,rollright,
        vytop0,vxbottom0
        
	vxleft_up      As for type_bc=3 family
	vxleft_down    As for type_bc=3 family
	yleft1         As for type_bc=3 family
	yleft2         As for type_bc=3 family
	rollleft       As for type_bc=3 family
	vxright_up     As for type_bc=3 family
	vxright_down   As for type_bc=3 family
	yright1        As for type_bc=3 family 
	yright2        As for type_bc=3 family
	rollright      As for type_bc=3 family
	vytop0         As for type_bc=3 family
	vxbottom0      As for type_bc=3 family      


                           **** CAUTION ****
      
        type_bc = 3 has special (likely inconsistent) properties when
        used with roller side boundary conditions.
        In the tectonic calculations all nodes on the side boundaries
        above the base are rollers, therefore the surface nodes on
        the sides can move vertically.
        
        However, later when the surface is updated for tectonic
        advection, the vertical advection of the surface nodes at
        the sides is set to 0, thereby suppressing their advection
        (inconsistent).
        
        The advection of the surface nodes at the sides due to
        surface processes is calcuated independently and depends on
        the parameter nerode_ends.
        
        Note that all other members of the type_bc=3 family allow
        vertical advection of the surface nodes at the sides. 
        type_bc=3 was the original version and should only be used in
        special circumstances.

        [2003/02/10]  Use type_bc = 300 instead of type_bc = 3 if you
        want to use this kind of boundary condition.  Type_bc = 300 is
        a corrected version of type_bc = 3.  Its input variables are
        exactly the same as for type_bc = 3.

           ----------------------------------------------        		      

          if(type_bc.eq.4)then
            read(8,*)vxleft_up,vxleft_down,yleft1,yleft2,rollleft,
                vxright_up,vxright_down,yright1,yright2,rollright,
                vytop0,vxbottom0,vxtop0

        Type_bc 4 is based on type_bc 3, with (1) corrected behaviour for
        side roller boundaries (the surface endpoints are allowed to move)
        and (2) the ability to restrict horizontal movement at the surface.
        This means that with type_bc=4, we can totally restrict movement
        at the surface (i.e. if vxtop0 = 1 and vytop0=1).

	vxleft_up      As for type_bc=3 family
	vxleft_down    As for type_bc=3 family
	yleft1         As for type_bc=3 family
	yleft2         As for type_bc=3 family
	rollleft       As for type_bc=3 family
	vxright_up     As for type_bc=3 family
	vxright_down   As for type_bc=3 family
	yright1        As for type_bc=3 family 
	yright2        As for type_bc=3 family
	rollright      As for type_bc=3 family
	vytop0         As for type_bc=3 family
	vxbottom0      As for type_bc=3 family      
        vxtop0    = 1  no horizontal velocities on top boundary
                  = 0  horizontal velocities allowed on top boundary


           ----------------------------------------------        		      

 

      if(type_bc.eq.400)then
      read(8,*)vpunchx,vpunchy,xsleft,xsright
      endif

      if(type_bc.eq.5)then
      read(8,*)ppunchx,ppunchy,xsleft,xsright
      endif

      if(type_bc.eq.6)then
      read(8,*)ppunchx,xpunch,ppunchy,ypunch,xsleft,xsright
      endif

      if(type_bc.eq.7)then
      read(8,*)ppunchx
      endif

      if(type_bc.eq.8)then
      read(8,*)omegarot,domega,xsleft,xsright
      endif


      IF type_bc = 10 then
      read(8,*)vleft,vright,xsleft,xsright,
               ampnoise,wavlnoise,decaylnoise

         This is the same as type_bc=1 with the addition of 
         sinusoidal noise to the x velocities between vleft and vright.
             
           vleft = velocity at far left of singularity
           vright = velocity at far right of singularity
           xsleft   }  [see diagram]  velocity changes linearly from vleft
           xsright  }                 to vright across the zone
           
               xsleft            xsright
                 |                 |
               --+--------o--------+---  
                 |        S        |


            ampnoise = amplitude of noise as a velocity in SI, m/s
            wavlnoise = wavelength of sinusoidal noise in SI, m
            decaylnoise = length scale in SI, m, over which the noise decays
                          to zero adjacent to xsleft and xsright
 
	 Choice of noise amplitude. For the noise velocity to give a
	 superimposed noise strain rate equal to alpha x mean velocity
	 gradient (where mean velocity = epsilon xx tectonic strain rate),
	 ampnoise should be set to:
 
             ampnoise = alpha x ( deltavx/ltrans)  x ( wavlnoise/4 )
 
         where deltavx = difference of velocity (m/s) across entire width of
                         transition zone (vxright-vxleft)
               ltrans = width (m) of transition zone (xsright-xsleft)
               wavelnoise = wavelength of the noise (m)

      IF type_bc = 11 then

           type_bc = 11 for singularity model   --------o---------
                                                        S

      read(8,*)vleft,rollleft, vright, rollright, xsleft,xsright
           vleft = velocity at far left of singularity
	   rollleft
	      1 = allow left endpoint of surface to move vertically
	      0 = don't allow left endpoint of surface to move vertically
           vright = velocity at far right of singularity
	   rollright
	      1 = allow right endpoint of surface to move vertically
	      0 = don't allow right endpoint of surface to move vertically
           xsleft   }  [see diagram]  velocity changes linearly from vleft
           xsright  }                 to vright across the zone
           
               xsleft            xsright
                 |                 |
               --+--------o--------+---  
                 |        S        |
          

                        WARNING! 
	   
	   type_bc=11 allows vertical rolling of the
           the left and right boundaries, except at the base.
           No modification is made to the Lagrangian grid
           outside the Eulerian domain.  This means there will be
           a mismatch between the Lagrangian and Eulerian meshes
           at boundaries where rolling has occurred.  This is
           probably okay for cases where material is advected out
           of the Eulerian domain.  It is also okay for boundaries
           where there is no advection into the domain.  It is
           up to the user to make sure this is satisfied.

 
      IF type_bc = 12 then
      
      This is like type_bc=11, but allows x velocities to vary along
      side boundaries (as type_bc 3 does).

      read(8,*) vxleft_up, vxleft_down, yleft1, yleft2, rollleft,
                vxright_up, vxright_down, yright1, yright2, rollright,
                xsleft,xsright

	   vxleft_up = x velocity at left boundary above yleft1
	   vxleft_down = x velocity at left boundary below yleft2

		 In the area between yleft1 and yleft2, the
		 boundary condition makes a smooth transition between
		 vxleft_up and vxleft_down.

	   yleft1 = vertical position at left boundary, above which
	            the x velocity is vxleft_up
	   yleft2 = vertical position at left boundary, below which
	            the x velocity is yleft_down
	   rollleft
	      1 = allow left endpoint of surface to move vertically
	      0 = don't allow left endpoint of surface to move vertically

	   vxright_up = x velocity at right boundary above yright1
	   vxright_down = x velocity at right boundary below yright2

		 In the area between yright1 and yright2, the
		 boundary condition makes a smooth transition between
		 vxright_up and vxright_down.

	   yright1 = vertical position at right boundary, above which
	            the x velocity is vxright_up
	   yright2 = vertical position at right boundary, below which
	            the x velocity is yright_down
	   rollright
	      1 = allow right endpoint of surface to move vertically
	      0 = don't allow right endpoint of surface to move vertically
	      
           xsleft   }  [see diagram]  velocity changes linearly from vleft
           xsright  }                 to vright across the zone
           
               xsleft            xsright
                 |                 |
               --+--------o--------+---  
                 |        S        |
          

                        WARNING! 
	   
	   type_bc=12 allows vertical rolling of the
           the left and right boundaries, except at the base.
           No modification is made to the Lagrangian grid
           outside the Eulerian domain.  This means there will be
           a mismatch between the Lagrangian and Eulerian meshes
           at boundaries where rolling has occurred.  This is
           probably okay for cases where material is advected out
           of the Eulerian domain.  It is also okay for boundaries
           where there is no advection into the domain.  It is
           up to the user to make sure this is satisfied.

 

      IF type_bc = 30 then
      read(8,*)vxleft_up,vxleft_down,yleft1,yleft2,rollleft,
                vxright_up,vxright_down,yright1,yright2,rollright,
                vytop0, vxbottom0,vybottom,decaylvybottom,
                ampnoise,wavlnoise,decaylnoise,xsleft,xsright

          type_bc = 30 is like type_bc=3, but with added sinusoidal noise
          in the x velocity along the base.

          NOTE that choosing this bc means you have a
          no-slip bottom boundary condition, no matter
          what vxbottom0 is set to.
          
	vxleft_up      As for type_bc=3 family
	vxleft_down    As for type_bc=3 family
	yleft1         As for type_bc=3 family
	yleft2         As for type_bc=3 family
	rollleft       As for type_bc=3 family
	vxright_up     As for type_bc=3 family
	vxright_down   As for type_bc=3 family
	yright1        As for type_bc=3 family 
	yright2        As for type_bc=3 family
	rollright      As for type_bc=3 family
	vytop0         As for type_bc=3 family
	vxbottom0      Value of vxbottom0 is ignored here, since
                       choosing this boundary condition means you have a
                       no-slip bottom boundary.
       	vybottom       As for type_bc=3 family
       	decaylvybottom As for type_bc=3 family
	xvyleft        As for type_bc=3 family
	xvyright       As for type_bc=3 family
	ampnoise = amplitude of noise as a velocity in SI, m/s
	wavlnoise = wavelength of sinusoidal noise in SI, m
	decaylnoise = length scale in SI, m, over which the noise decays
                      to zero adjacent to xsleft and xsright
	xsleft= left limit in SI, m, for the added basal noise function
	        and for y velocity on bottom boundary
	xsright= right limit in SI, m, for the added basal noise function
	        and for y velocity on bottom boundary


      if(type_bc.eq.31)then
        read(8,*)vxleft_up,vxleft_down,yleft1,yleft2,rollleft,
     &           vxright_up,vxright_down,yright1,yright2,rollright,
     &           vytop0, vxbottom0,vybottom,decaylvybottom,
     &           ampnoise,decaylnoise,xsleft,xsright

          type_bc = 31 is like type_bc=3, but with added white noise
          in the x velocity along the base.
 
          NOTE that choosing this bc means you have a
          no-slip bottom boundary condition, no matter
          what vxbottom0 is set to.

	  vxleft_up      As for type_bc=3 family
	  vxleft_down    As for type_bc=3 family
	  yleft1         As for type_bc=3 family
	  yleft2         As for type_bc=3 family
	  rollleft       As for type_bc=3 family
	  vxright_up     As for type_bc=3 family
	  vxright_down   As for type_bc=3 family
	  yright1        As for type_bc=3 family 
	  yright2        As for type_bc=3 family
	  rollright      As for type_bc=3 family
	  vytop0         As for type_bc=3 family
	  vxbottom0      Value of vxbottom0 is ignored here, since
                	 choosing this boundary condition means you have a
                	 no-slip bottom boundary.
	  vybottom       As for type_bc=3 family
       	  decaylvybottom As for type_bc=3 family
	  xvyleft        As for type_bc=3 family
	  xvyright       As for type_bc=3 family
	  ampnoise = amplitude of noise as a velocity in SI, m/s
	  decaylnoise = length scale in SI, m, over which the noise decays
                	to zero adjacent to xsleft and xsright
	  xsleft= left limit in SI, m, for the added basal noise function
	        and for y velocity on bottom boundary
	  xsright= right limit in SI, m, for the added basal noise function
	        and for y velocity on bottom boundary


      IF type_bc = 312 then
      read(8,*)vxleft_up,vxleft_down,yleft1,yleft2,rollleft,
               vxright_up,vxright_down,yright1,yright2,rollright,
               vytop0, vxbottom0,vybottom,decaylvybottom,
               ampnoise,decaylnoise,xsleft,xsright
     
          type_bc = 312 is like type_bc = 3, but with added white noise
          in the y velocity along the base.

          NOTE that choosing this bc means you can have a
          free-slip bottom boundary condition - an improvement
          from type_bc.eq.31!

	  vxleft_up      As for type_bc=3 family
	  vxleft_down    As for type_bc=3 family
	  yleft1         As for type_bc=3 family
	  yleft2         As for type_bc=3 family
	  rollleft       As for type_bc=3 family
	  vxright_up     As for type_bc=3 family
	  vxright_down   As for type_bc=3 family
	  yright1        As for type_bc=3 family 
	  yright2        As for type_bc=3 family
	  rollright      As for type_bc=3 family
	  vytop0         As for type_bc=3 family
	  vxbottom0      As for type_bc=3 family
	  vybottom       As for type_bc=3 family
       	  decaylvybottom As for type_bc=3 family
	  xvyleft        As for type_bc=3 family
	  xvyright       As for type_bc=3 family
	  ampnoise = amplitude of noise as a velocity in SI, m/s
	  decaylnoise = length scale in SI, m, over which the noise decays
                	to zero adjacent to xsleft and xsright
	  xsleft= left limit in SI, m, for the added basal noise function
	        and for y velocity on bottom boundary
	  xsright= right limit in SI, m, for the added basal noise function
	        and for y velocity on bottom boundary


      IF type_bc = 32 then
      read(8,*)vxleft_up,vxleft_down,yleft1,yleft2,rollleft,
               vxright_up,vxright_down,yright1,yright2,rollright,
               vytop0, vxbottom0,vybottom,decaylvybottom,
               ampnoise,wavlnoise,decaylnoise,xsleft,xsright,
               noisenodey1,noisenodey2
     
	  type_bc = 32 is like type_bc = 3, but with added sinusoidal
	  noise in the x body force along the nodes in rows noisenodey1
	  to noisenodey2 (inclusive).

	  vxleft_up      As for type_bc=3 family
	  vxleft_down    As for type_bc=3 family
	  yleft1         As for type_bc=3 family
	  yleft2         As for type_bc=3 family
	  rollleft       As for type_bc=3 family
	  vxright_up     As for type_bc=3 family
	  vxright_down   As for type_bc=3 family
	  yright1        As for type_bc=3 family
	  yright2        As for type_bc=3 family
	  rollright      As for type_bc=3 family
	  vytop0         As for type_bc=3 family
	  vxbottom0      As for type_bc=3 family
	  vybottom       As for type_bc=3 family
          decaylvybottom As for type_bc=3 family
	  xvyleft        As for type_bc=3 family
	  xvyright       As for type_bc=3 family
          ampnoise = amplitude of noise as a velocity in SI, m/s
          wavlnoise = wavelength of sinusoidal noise in SI, m
          decaylnoise = length scale in SI, m, over which the noise decays
                        to zero adjacent to xsleft and xsright
          xsleft= left limit in SI, m, for the added noise function on
                  rows noisenodey1 to noisenodey2
                   and for y velocity on bottom boundary
          xsright= right limit in SI, m, for the added noise function on
                   rows noisenodey1 to noisenodey2
                   and for y velocity on bottom boundary
          noisenodey1 = top row number where noise is to be added to x body 
                        force
          noisenodey2 = bottom row number where noise is to be added to x body 
                        force



      IF type_bc = 33 then
      read(8,*)vxleft_up,vxleft_down,yleft1,yleft2,rollleft,
                vxright_up,vxright_down,yright1,yright2,rollright,
                vytop0, vxbottom0, vybottom,decaylvybottom,
                ampnoise,decaylnoise,xsleft,xsright,
                noisenodey1, noisenodey2
     
	  type_bc = 33 is like type_bc = 3, but with added white
	  noise in the x body force along the nodes in rows noisenodey1
	  to noisenodey2 (inclusive).

	  vxleft_up      As for type_bc=3 family
	  vxleft_down    As for type_bc=3 family
	  yleft1         As for type_bc=3 family
	  yleft2         As for type_bc=3 family
	  rollleft       As for type_bc=3 family
	  vxright_up     As for type_bc=3 family
	  vxright_down   As for type_bc=3 family
	  yright1        As for type_bc=3 family
	  yright2        As for type_bc=3 family
	  rollright      As for type_bc=3 family
	  vytop0         As for type_bc=3 family
	  vxbottom0      As for type_bc=3 family
	  vybottom       As for type_bc=3 family
          decaylvybottom As for type_bc=3 family
	  xvyleft        As for type_bc=3 family
	  xvyright       As for type_bc=3 family
          ampnoise = amplitude of noise as a velocity in SI, m/s
          decaylnoise = length scale in SI, m, over which the noise decays
                        to zero adjacent to xsleft and xsright
          xsleft = left limit in SI, m, for the added noise function on
                  rows noisenodey1 to noisenodey2
                   and for y velocity on bottom boundary
          xsright = right limit in SI, m, for the added noise function on
                  rows noisenodey1 to noisenodey2
                   and for y velocity on bottom boundary
          noisenodey1 = top row number where noise is to be added to x 
                        body force
          noisenodey2 = bottom row number where noise is to be added to x 
                        body force


      [Eulerian remeshing]
      read(8,*)ksplitx1,tsmooth,ismooth, leftsmooth, rightsmooth
            ksplitx1 = number of subintervals per Eulerian delta X
                       at which markers are embedded for top surface
                       tracking (use 2, or < 4)
            tsmooth = number of timesteps between two numerical smoothings             
		      of the top surface, i.e. smoothing is done
		      every tsmooth timesteps.  NOTE: in the timestep
		      when smoothing is done, smoothing is done ncyc2
		      times (see subroutine vylifts).

		                     ** CAUTION ***
		                     
		      Surface smoothing may not be consistent with
		      all types of tectonic and surface processes. In
		      earlier versions surface smoothing is applied
		      to the vertically advected Eulerian tectonic
		      surface and then surface and other processes
		      are applied to the smoothed surface.
		      
		      in the nov21-02 version of the code, surface
		      smoothing is applied to the final surface after
		      the effects of advection owing to the surface
		      processes have been added.
		      
		      Numerical surface smoothing is best avoided if
		      possible.  Instead, it is preferable to use the
		      diffusion option in the surface processes - see
		      subroutine erode.  If smoothing is all that is
		      required (not true surface processes) ensure
		      that the material properties of any injected
		      Lagrangian particles that correspond to
		      sediment are the same as those of the material
		      that comprises the surface that is being
		      smoothed.

		      
            ismooth = number of smoothing passes per smoothing timestep
                      (= 1)
                      
            leftsmooth, rightsmooth = left and right columns which define
		      the portion of the surface to be smoothed.  To smooth
		      the entire surface (the usual procedure), set both
		      leftsmooth and rightsmooth to 0 (or, you could set
		      leftsmooth to 2 and rightsmooth to nx-1).
      
      read(8,*)vscale,verror,verror1,miterv,irescale,kinitvis,
	       force_iter1, force_iter, force_continue,
	       mech_trouble_iter 

            vscale = velocity scale for scaling velocity errors in
                     Picard iterations.  Convergence is defined to be
		     achieved when 
			 maxdv/vscale .le. verror
                     where maxdv is the greatest change in velocity
		     (in either x or y direction) since the previous
		     iteration.

            verror = relative error for timesteps > 1 in Picard iteration
            verror1 = relative error for timestep 1 in Picard iterations
            miterv = max number of iterations in Picard iterations. If
		     convergence is not achieved by miterv iterations,
		     the model will stop unless force_continue is set
		     greater than 0.
            irescale = scheme used for updating viscosities
                       0 = calculate viscosities according to material
                           properties if stress above yield viscosities
                           are rescaled
                           -> viscosity rescaling created by strain rate
                           (see earlier notes)
            kinitvis = max number of iterations per step before the
                       viscosity is reinitialized and computation is started
                       from scratch.  Should always check by setting
                       kinitvis = 1.
            
	    force_iter1 = number of iterations to force for the first
			  timestep.  This number of iterations will
			  always be performed even if convergence is
			  achieved earlier.  If convergence has not been
			  reached after the forced iterations, the code
			  will continue to iterate until convergence is
			  reached or until the maximum iterations (miterv)
			  is reached.
			  0 = don't force a minumum number of iterations
			      (normal behaviour)
	    force_iter = number of iterations to force for timesteps > 1
			 in the Picard iterations.  This number of
			 iterations will always be performed even if
			 convergence is achieved earlier.  If convergence
			 has not been reached after the forced iterations,
			 the code will continue to iterate until
			 convergence is reached or until the maximum
			 iterations (miterv) is reached.
			 0 = don't force a minimum number of iterations.
			     (normal behaviour)
	    force_continue = flag to indication whether to force continued
			     timesteps when Picard divergence occurs.
			     0 = if Picard divergence occurs, stop the model
				 (normal behaviour)
			     1 = if Picard divergence occurs, don't stop
				 the model but continue to next timestep
				 anyway.  USE WITH CAUTION! ONLY TO BE USED
				 IN SPECIAL CASES!
            
	    mech_trouble_iter = iteration number at which to write an output
				file.  This output file is given frame number
				99, which means that its name will be something
				like "g01_p00_f99_o" or "g02_p00_f99_o". This
				file is useful for debugging problem models
				where convergence is not reached on the last
				timestep and the restart file doesn't give any
				clues why.
				If mech_trouble_iter is not specified (as in
				earlier versions of SOPALE1_i), it will be set
				to a default value of 50 iterations.
          
      
      read(8,*)ncyc2,color2op,lagtest, eulcolorrule
            ncyc2 = number of cycles in Lag/Eul remeshing/tracking per step
                    (use 2)
                    includes top surface recycling (0 gives no tracking)
            color2op = 0 -> Lagrangian particles color the Eulerian grid
                            when computation initialized
                            (Microfem-like process)
                       1 -> Eulerian grid colors Lagrangian particles
                            initially
                       2 -> Philippe's code says this means
                            "nonmixedelementcoloring grid 1 with
                            boxes/strings"                            
		       10  -> read intial Lagrangian particle colors
		              (materials) from file lagcolors_i
		              
            lagtest - use 0
                      1 -> tracking with uniform velocity field, useful
                           for benchmarks
            
            eulcolorrule = rule to use for coloring Eulerian cells from
                           Lagrangian particles.
                           0 => use old method: color Eulerian cell the
                                same as the last Lagrangian particle found
                                in the cell.  Note that this will bias the
                                material colors in favour of Lagrangian
                                particles with higher node number.
                           1 => use "majority rules": Eulerian cell is colored
                                according to the majority (actually plurality)
                                of Lagrangian particles; i.e. the color which
                                is represented by the most number of Lagrangian
                                particles in the Eulerian cell.  If there is
                                not a majority, then the Eulerian cell's 
                                existing color is not changed.
                 
                       
      read(8,*)ucompute,tcompute,ale_compute,dt_compute,benchmark,
               tbottom,ttop,e_temperature,cfl,
	       lambda_blend1, lambda_blend2 
      
            ucompute = 1 -> compute a velocity using an FE computation
                     = 0 -> no FE computation, computation of another sort
            tcompute = 1 -> compute a temperature field using an FE 
                            computation
                     = 0 -> no FE computation
                            0 for no thermo-mech
            ale_compute = 1 -> ALE computation with remeshing and computation
                        = 0 -> no ALE computation /no remeshing
            dtcompute   = 0 -> not going to recompute dt
                        = 1 -> CFL recomputation of dt
            tbottom = bottom temperature (Absolute temperature, not Celsius)
                      Linear geotherm top to bottom, no radioactive heating
            ttop = top temperature (absolute temperature)
            e_temperature = 0 -> Eulerian temps that control material props
                                 are elemental (one value only)
                          = 1 -> Eulerian temps that control material props
                                 are nodal
                          Should use 0 when temperatures are not computed,
                          i.e. they are just advected by markers, i.e.
                           Lag node ----------------> Eulerian element
                             |                          ^
                             |                          |
                             v                          |
                            Average of lag nodes => ----|
                           
                           When no lag nodes in element, keep value
                           Temperatures are not updated
                           
                           Entity is 'per cell' level, i.e. each element has
                           only one property and properties are not specified
                           at the Gauss point level
            cfl = used to automatically adjust timestep length (dt) when
                  dt_compute = 1.
                  Timestep length is reset so that the cfl takes value given
                  in cfl.


	    lambda_blend1, lambda_blend2
		= interval during the first timestep where you want to 
		  "blend in" the value of lambda_pore_prs (fluid pore 
		  pressure ratio) in the calculation of plastic yielding 
		  (splas).
		  In timestep 1, starting at iteration number lambda_blend1,
		  the value of lambda_pore_prs is gradually increased from
		  zero until it reaches its full value at iteration number
		  lambda_blend2.
		  This technique avoids problems with negative values of
		  plastic yielding in the first timestep, when the model
		  has not yet stabilized.  Normally used in conjunction
		  with plas_fluid_method (see above).
		  NOTE: You should make sure that force_iter1 is at least
		  as large as lambda_blend2.



      
      read(8,*)steady,analy,esteady
            steady = 1 for steady state computation
                   = 0 for non-steady-state computation
                       (steady state implies dt can be 0)
            analy = obsolete, set to 0
            esteady = only used when steady = 1
                    = the error level you want to wait for until you reach
                      steady state
      
      *** Start thermal section (only needed if tcompute = 1) ****
      
      read(8,*)maxter,dtempm,ther_advection,th_gamma_time,th_cn,tempsmooth
            maxter  = maximum error for thermal iterations
		      For steady=0, this is a relative error
		      For steady=2, this is an absolute error
            dtempm  = temperature error (units are temperature). If temperature
                      error more than this, continue iterating.
		      MAY NOT BE USED (value is set in code)
            ther_advection = 1 there is advection of the temperature field
                           = 0 there is no advection of the temperature field
            th_gamma_time  = numerical parameter for integration in time
                             Choose between .5 and 1, typically 1
                             Trapezoidal schemes for time integration
            th_cn     =  Krank-Nicholson parameter
                         Use 1
            tempsmooth = timestep interval between successive smoothings of
			 temperature field
      
      read(8,*)type_bct
            type_bct   = boundary condition type for thermal problem
                         DURING the model run (not initialization).
			 See type_bct0 for initial conditions in dmtm.
                         1 = convection
                         2 = convection
                         3 = for subduction models
			     for deep model thermo-mechanical (dmtm)
                         
         If type_bct = 3, next line defines bct3 parameters (11 of them).
	 This line is not needed if type_bct is not 3.
           bct3(1) = top temperature (degrees K) [see tshift in the code]
           bct3(2) = 0 bottom is imposed temperature, value bct3(3)
                   = 1 bottom is imposed flux, value bct3(3)
           bct3(3) = temperature or flux at base (depending on value of
                      bct3(2)
           bct3(4) = 0  On left boundary, impose initial geotherm above
                        bct3(5) and bct3(7) for flux below bct3(5)
                   = 1  On left boundary, impose bct3(6) as temperature
                        above bct3(5) and bct3(7) as flux below bct3(5)
                   = 2  On left boundary, impose bct3(6) as temperature
                        above bct3(5) and initial geotherm below bct3(5)
           bct3(5) = z  Position on left boundary about which the boundary 
			condition changes (see bct3(4)) - distance in metres
			or user units measured from origin upwards
           bct3(6) = T  Temperature on left boundary when bct3(4) =1 or 2
			Not used when bct3(4) = 0.
           bct3(7) = flux  Lateral flux imposed at left boundary below
			   bct3(5) when bct3(4) = 0 or 1.   Not used when
			   bct3(4) = 2.
           bct3(8) = 0  On right boundary, impose initial geotherm above
                        bct3(9) and bct3(11) for flux below bct3(9)
                   = 1  On right boundary, impose bct3(10) as temperature
                        above bct3(9) and bct3(11) as flux below bct3(9)
                   = 2  On right boundary, impose bct3(10) as temperature
                        above bct3(9) and initial geotherm below bct3(9)
           bct3(9) = z  Position on right boundary about which the boundary 
			condition changes (see bct3(8)) - distance in metres
			or user units measured from origin upwards
           bct3(10) = T Temperature on right boundary when bct3(8) =1 or 2
			Not used when bct3(8) = 0.
           bct3(11) = flux  Lateral flux imposed at right boundary below
			    bct3(9) when bct3(8) = 0 or 1.   Not used when
			    bct3(8) = 2.
      
      
      read(8,*)th_dt,cfl,vzgiven
            th_dt = timestep for thermal code
		    Is computed internally when steady=0 and dtcompute=1
                    (steady=0 means non steady state calculation)
            cfl
                     NOTE: cfl is re-read; the value used is the last value
                     taken
            vzgiven = only used for benchmarks 20, 21, 22, 23
		      not used in deep thermo-mechanical models
            
      read(8,*)ncolor2t, added_sed_colort
            ncolor2t  = number of thermal materials
	    added_sed_colort = thermal material number
	                       to use for injected sediment particles

      For each thermal material (1 to ncolor2t), have the following 5 lines:
         read(8,*)thermod1_conductivity(i)
                = K, thermal conductivity
		units W/(m K) or user units
      
         read(8,*)thermod1_rau0(i)
                = Thermal density
                  Density used by thermal code
		  units kg/m**3 or user units
      
         read(8,*)thermod1_cp(i)
                = Specific heat
		  units J/(kg K) or user units
      
         read(8,*)thermod1_fheating(i)
		= frictional heating coefficient
                 Not active, use 0
      
         read(8,*)thermod1_power_permass(i)
		= internal heat production (e.g. radiactive heat generation)
                 units watts/kg
                (NOTE: need to convert from A in units of watts/m**3 to
		       H in units of watts/kg.  A/rau0 = H)


      *** End thermal section (only needed if tcompute = 1) ****
      
      READ(8,*) OUTFIL
            OUTFIL = name of output file in which blkfct will write 
                     diagnostic info
                     
            Per linear SPD system (i.e. per blkfct system) need 4 lines
            for each of MECH and THERMAL SYSTEM if they are needed
            
      READ(8,*) MATFIL
            MATFIL = file name containing profile of the matrix to
                     be solved
                     
                     Can have 0, 1, or 2 systems to solve
                     
                     These are all in "matrices subdirectory" and
                     are generated by program pmain.f in subdirectory
                     "sources".

     READ(8,*) ICASE
           ICASE = 1 natural ordering
                 = 2 mmd ordering
                 
     READ(8,*) CACHSZ
           CACHSZ = cachesize in kilobytes
                    Try 64, can be flexible
                    
     READ(8,*) LEVEL
           LEVEL = loop unrolling level
                   Try 8

     *** Begin thermal section (only needed if tcompute = 1) **** 

            
      READ(8,*) MATFIL (for thermal systems)
            MATFIL = file name containing profile of the matrix to
                     be solved
                     
                     Can have 0, 1, or 2 systems to solve
                     
                     These are all in "matrices subdirectory" and
                     are generated by program pmain.f in subdirectory
                     "sources".

     READ(8,*) ICASE (for thermal systems)
           ICASE = 1 natural ordering
                 = 2 mmd ordering
                 
     READ(8,*) CACHSZ (for thermal systems)
           CACHSZ = cachesize in kilobytes
                    Try 64, can be flexible
                    
     READ(8,*) LEVEL (for thermal systems)
           LEVEL = loop unrolling level
                   Try 8
 
     *** End thermal section (only needed if tcompute = 1) ***
      
      read(8,*)ibinread
            ibinread = 1   binary read/writes to outputs
                     = 0   formatted read/writes (may not work)
      
      NOTE: the following filenames are compulsory in terms of number
      of characters, i.e. the filenames cannot have any more or less
      characters than the format specified in the read statement.
      Also, the first 13 characters of the filename are assumed to 
      specify the directory (including the trailing "/").
      
      read(8,'(a26)')filename1
            filename1 = name of restart file for Eulerian grid
 
      read(8,'(a26)')filename2
            filename2 = name of restart file for Lagrangian grid

      read(8,'(a33)')filenameg1
            filenameg1 = name of file to save timestep information
                         
      read(8,'(a26)')filenameg3
            filenameg3 = name of file containing grid type and topology
                         used by mozart_plot.x for Eulerian 

      read(8,'(a26)')filenameg4
            filenameg4 = name of file containing grid type and topology
                         used by mozart_plot.x for Lagrangian 
 
  *** Mechanical material (colour) boxes section ****

      NOTE: THIS INPUT DEPENDS ON color2op.  If color2op = 0, we will 
      read these lines
                    
      read(8,*)ncolor2
            ncolor2 = number of Lagrangian mechanical material boxes 
		      defined on lag grid
            
            Then for each box, read the 4 vertices (in metres or user units)
	    of the box in counterclockwise order, i.e.:

            read(8,*)xx1,xx2,xx3,xx4
                  xx1 =  x coordinate of top left corner of box
                  xx2 =  z coordinate of top left corner of box
                  xx3 =  x coordinate of bottom left corner
                  xx4 =  z coordinate of bottom left corner
            read(8,*)xx5,xx6,xx7,xx8
                  xx5 =  x coordinate of bottom right corner
                  xx6 =  z coordinate of bottom right corner
                  xx7 =  x coordinate of top right corner
                  xx8 =  z coordinate of top right corner 
            
            read(8,*)newcol2
                  newcol2 = corresponding material number for material
			    particles in this box
      
      If color2op = 1 then define the layers on Eulerian grid :

          read(8,*)nelayers
                nelayers = number of layers
                
                For each layer define the last row of the Eulerian grid
		for this layer, counting top to bottom
                
                read(8,*)(layerse(i),i=1,nelayers)
                    layerse(i) are the last row on the Eulerian grid
		    for each of the nelayers
                      
                Each layer is automatically given the mechanical material 
		(colour) number in order of the layers, i.e. layer 1 is
		given material number 1, layer 2 is material 2, etc.
                Round off error determines box.
                

      NOTE: Slight perturbation on lag grid = 0.333 x 10e-6 in x
      

      Mixing model determines material property number to be taken
      by an Eulerian element during regridding when there are Lagrangian
      particles with different material property numbers in the Eulerian
      element. The mixing model either uses the minimum material property
      number (color number), or the material property number of the last
      particle found (*** needs to be checked ***).
                    
   
  *** Thermal material (colour) boxes section *** 
    
      DMTM version note: You need to define at least one thermal material
      box, even if you're not doing thermal calculations.  If you are
      doing a mechanical-only model, define a dummy thermal material
      box which covers your entire domain and consists of thermal
      material 1. (You don't have to define a thermal material).      

                    
      read(8,*)ncolor2
            ncolor2 = number of Lagrangian thermal material boxes 
		      defined on lag grid
            
            Then for each box, read the 4 vertices (in metres or user units)
	    of the box in counterclockwise order, i.e.:

            read(8,*)xx1,xx2,xx3,xx4
                  xx1 =  x coordinate of top left corner of box
                  xx2 =  z coordinate of top left corner of box
                  xx3 =  x coordinate of bottom left corner
                  xx4 =  z coordinate of bottom left corner
            read(8,*)xx5,xx6,xx7,xx8
                  xx5 =  x coordinate of bottom right corner
                  xx6 =  z coordinate of bottom right corner
                  xx7 =  x coordinate of top right corner
                  xx8 =  z coordinate of top right corner 
            
            read(8,*)newcol2
                  newcol2 = corresponding thermal material number for material
			    particles in this box
      
  *** End of thermal material (colour) boxes section *** 



      read(8,*)dmtm, n_skip_thermal
            dmtm = 1  flag for deep thermo-mechanical model
            dmtm = 0  if you are doing mechanical model only
            n_skip_thermal = number of iterations where the thermal
                             calculations are skipped (i.e. if
                             n_skip_thermal = 50, then the code will
                             not perform thermal calculations for the 
                             first 50 iterations)
      
      Only if dmtm = 1 and rres is NOT 1, then read the following initial
      thermal boundary conditions: 

      read(8,*)type_bct0

            type_bct0   = initial thermal boundary condition type
			  for deep model thermo-mechanical (dmtm)

         If type_bct0 = 0, then the temperature field as defined in
                	   sopale1_plus_i is used as the initial
			   temperature field. (In this case, the code
			   skips over the finite element steady-state
			   solution.)

         If type_bct0 = 3, next line defines bct30 parameters (11 of them).
			   The code sets these initial boundary conditions
			   and does a finite element steady-state solution.

           read(8,*)(bct30(i),i=1,11)
 	 
           bct30(1) = top temperature (degrees K) [see tshift for the code]
           bct30(2) = 0 bottom is imposed temperature, value bct30(3)
                    = 1 bottom is imposed flux, value bct30(3)
           bct30(3) = temperature or flux at base (depending on value of
                      bct30(2)
           bct30(4) = 0  On left boundary, impose initial geotherm above
                         bct30(5) and bct30(7) for flux below bct30(5)
                    = 1  On left boundary, impose bct30(6) as temperature
                         above bct30(5) and bct30(7) as flux below
                         bct30(5)
                    = 2  On left boundary, impose bct30(6) as temperature
                         above bct30(5) and initial geotherm below
                         bct30(5)
           bct30(5) = z  Position on left boundary about which the boundary 
			 condition changes (see bct30(4)) - distance in metres
			 or user units measured from origin upwards
           bct30(6) = T  Temperature on left boundary when bct30(4) =1 or 2
			 Not used when bct30(4) = 0.
           bct30(7) = flux  Lateral flux imposed at left boundary below
			    bct30(5) when bct30(4) = 0 or 1.   Not used when
			    bct30(4) = 2.
           bct30(8) = 0  On right boundary, impose initial geotherm above
                         bct30(9) and bct30(11) for flux below bct30(9)
                    = 1  On right boundary, impose bct3(10) as temperature
                         above bct30(9) and bct30(11) as flux below
                         bct30(9)
                    = 2  On right boundary, impose bct30(10) as temperature
                         above bct30(9) and initial geotherm below
                         bct30(9)
           bct30(9) = z  Position on right boundary about which the boundary 
			 condition changes (see bct30(8)) - distance in metres
			 or user units measured from origin upwards
           bct30(10) = T Temperature on right boundary when bct30(8) =1 or 2
			 Not used when bct30(8) = 0.
           bct30(11) = flux  Lateral flux imposed at right boundary below
			    bct30(9) when bct30(8) = 0 or 1.   Not used when
			    bct30(8) = 2.

  
     
     NOTE no geotherm
     Full calculation of diffusion/production thermal equilibrium
     Note need for consistent temperature input - use Kelvin
  
      read(8,*)ntperts
            ntperts = number of boxes for temperature perturbation
      
      For each box, read in coordinates of corners and temperature
      perturbation value (tpert).  The value of tpert is added to
      the initial temperature values.

      do k=1,ntperts
        read(8,*)xx1,xx2,xx3,xx4
        read(8,*)xx5,xx6,xx7,xx8
        read(8,*)tpert

===========================================================================

File: sopale1_plus_i


This file contains parameters to define the Eulerian grid, the
properties of the temperature field on the Eulerian grid, the Lagrangian
grid, and the properties of the temperature field on the Lagrangian
grid.  It also defines any initial perturbation of "internal interfaces"
on the Eulerian and Lagrangian grids.  It's organized into a number of
sections; each section begins with an identifier.  The sections can be
in any order and there can be blank lines in between each section.

                          ** NOTE **
         
         (0,0) is defined to be the origin of the Eulerian grid.
         This means that Eulerian nodes will have only positive
         coordinates.  The Lagrangian grid, however, may have
         negative coordinates.
         
EGRID section
-------------

Defines the Eulerian grid. The first line of this section is the word
EGRID (on a line by itself).


      read(78,*)xor,lhoriz,nlayers,nbzeremesh
           xor      = x origin, in user units  [MUST BE ZERO]
           lhoriz   = grid width, in user units
           nlayers  = number of layers with different spacing in grid
	   nbzeremesh = row number in Eulerian grid, above which grid
			is uniformly dilated (stretched) during
			Eulerian regridding to account for changes in
			the top or bottom surfaces.  Generally chosen
			to match a change in the vertical refinement
			of the Eulerian grid. See subroutine
			zeremesh. When the whole grid is vertically
			dilated nbzeremesh = ny (total number of
			nodes vertically in the Eulerian grid.)
			Must be greater than the row number of the
			lowest interface (see row_interface below).

      read(78,*)(yl(i),i=1,nlayers+1)
           yl(i)    = depth (in user units) of the layer interfaces
                      (i.e. depths of the boundaries between the
                      layers).  These are defined from top to bottom,
                      where the bottom is at depth 0.

      read(78,*)(nyl(i),i=1,nlayers)
           nyl(i)    = number of elements in layer i
                       Note that the layers are numbered from top to
                       bottom.
      read(78,*) refh
           refh = reference height for erosion

      read(78,*) n_interface
           n_interface = number of interfaces (internal surfaces) which 
                         will be advected.  You must specify at least
                         one interface, which is the surface.
      do i=1, n_interface
          read(78,*) row_interface(i), apply_surfproc(i)
          
            row_interface(i) = row number of this interface. Rows are 
                               numbered from the top down.
            apply_surfproc(i) = flag to indicate whether to apply
                                surface processes (erosion, progradation,
                                etc.) to this surface
                                 0 => don't apply surface processes
                                 1 => apply surface processes


ETEMP section
-------------

Defines the initial Eulerian thermal structure. The first line of this section is
ETEMP (on a line by itself).

      read(78,*) nlt
           nlt = number of layers in the thermal grid
      
      read(78,*) ytop(i),ylthick(i),rada(i),cond(i),temptopK(i),qbase(i)
        ytop = position of top interface of layer in model coordinate system
        ylthick = thickness of layer
        rada = radioactive heat generation, A, ie production per unit volume
        cond = thermal conductivity of layer
        temptopK = temperature of upper interface of the layer (Kelvin)
        qbase = vertical heat flux into the base of the layer


LGRID section
-------------

Defines the Lagrangian grid. The first line of this section is the word
LGRID (on a line by itself)

      read(78,*)xor,lhoriz,nlayers
           xor      = x origin, in user units
           lhoriz   = grid width, in user units
           nlayers  = number of layers with different spacing in grid


      read(78,*)(yl(i),i=1,nlayers+1)
           yl(i)    = depth (in user units) of the layer interfaces
                      (i.e. depths of the boundaries between the
                      layers).  These are defined from top to bottom,
                      where the bottom is at depth 0.

      read(78,*)(nyl(i),i=1,nlayers)
           nyl(i)    = number of elements in layer i
                       Note that the layers are numbered from top to
                       bottom.

LTEMP section
-------------

Defines the initial Lagrangian thermal structure. The first line of this section
is the word LTEMP (on a line by itself).

      read(78,*) nlt
           nlt = number of layers in the thermal grid
      
      read(78,*) ytop(i),ylthick(i),rada(i),cond(i),temptopK(i),qbase(i)
        ytop = position of top interface of layer in model coordinate system
        ylthick = thickness of layer
        rada = radioactive heat generation, A, ie production per unit volume
        cond = thermal conductivity of layer
        temptopK = temperature of upper interface of the layer (Kelvin)
        qbase = vertical heat flux into the base of the layer

EPERTURB section
----------------

Defines the perturbation for the Eulerian grid, including any definition
of the initial surface.  If you don't want any perturbation for the
Eulerian grid, just leave this section out entirely.  (This section and
the LPERTURB section replace the old initsurf_i file.)

     read(10,*) nperturb
        nperturb = number of rows to perturb for the initial grid setup     
 
     Next lines:
     
      do i = 1, nperturb
        read(10,*) row_perturb(i), type_perturb(i)

          row_perturb(i) = row number (in grid) to perturb
          type_perturb(i) = type of perturbation to apply
 
          if type_perturb(i) = 1 (sinusoidal perturbation)
           read(10,*) xsleft_surf, xsright_surf, amp_surf, wavl_surf, 
                        decayl_surf,ymid_sin

           xsleft_surf  - left endpoint of sinuosoidal surface topography
           xsright_surf - right endpoint of sinuosoidal surface topography
           amp_surf     - amplitude of sine wave
           wavl_surf    - wavelength of sine wave
           decayl_surf  - length scale over which the amplitude of the
                          sine wave decays to zero adjacent to 
                          xsleft_surf and xsright_surf
           ymid_sin     - y position of midline of sine wave.  Choose
                          this value so that ymid_sin + amp_surf is
                          equal to or less than the height of the
                          Eulerian grid.
          
          if type_perturb(i) = 2 (Lehner Normal perturbation)
           read(10,*) h0, hinf, x0, prograde_length, h0star
	     h0         = surface height at x=0 (metres)
	     hinf       = limit of surface height at x = infinity (metres)
	     x0         = initial x position where surface height is at
	                   (h0-hinf)/2 + hinf (metres)
	     prograde_length = width of the distribution (metres) defined by
                               h(x) = (h0-hinf)exp(-((x-x0)**2)/prograde_length**2)
			        + hinf
	     h0star          = lowered height behind the prograding surface

          if type_perturb(i) = 3 (level-slope-level)
           read(10,*) x1, level1, x2, level2
             x1         = x position of start of slope section
	     level1     = y position of start of slope section
	     x2         = x position of end of slope section
	     level2     = y position of end of slope section
	     
	     For x < x1, y will be set to level1 
	     For x > x2, y will be set to level2
	     A straight-line segment joins (x1,level1) and (x2,level2).


     Next line:
     
     read(10,*)remeshbottom
        remeshbottom = row number in grid (counting from the top),
                     above which the vertical grid spacing is
                     dilated or compressed to accomodate the
                     perturbations specified; at this row and below,
                     the grid spacing is not affected.
                     
                     NOTE: make sure you choose remeshbottom 
                     to be further down than the lowest point
                     of your defined surface topography.


LPERTURB section
----------------

Defines the perturbation for the Lagrangian grid, including any definition
of the initial surface.  If you don't want any perturbation for the
Lagrangian grid, just leave this section out entirely.

     read(10,*) nperturb
        nperturb = number of rows to perturb for the initial grid setup     
 
     Next lines:
     
      do i = 1, nperturb
        read(10,*) row_perturb(i), type_perturb(i)

          row_perturb(i) = row number (in grid) to perturb
          type_perturb(i) = type of perturbation to apply
 
          if type_perturb(i) = 1 (sinusoidal perturbation)
           read(10,*) xsleft_surf, xsright_surf, amp_surf, wavl_surf, 
                        decayl_surf,ymid_sin

           xsleft_surf  - left endpoint of sinuosoidal surface topography
           xsright_surf - right endpoint of sinuosoidal surface topography
           amp_surf     - amplitude of sine wave
           wavl_surf    - wavelength of sine wave
           decayl_surf  - length scale over which the amplitude of the
                          sine wave decays to zero adjacent to 
                          xsleft_surf and xsright_surf
           ymid_sin     - y position of midline of sine wave.  Choose
                          this value so that ymid_sin + amp_surf is
                          equal to or less than the height of the
                          Eulerian grid.
          
          if type_perturb(i) = 2 (Normal perturbation)
           read(10,*) h0, hinf, x0, prograde_length, h0star
          
          if type_perturb(i) = 3 (level-slope-level)
           read(10,*) x1, level1, x2, level2

     Next line:
     
     read(10,*)remeshbottom
        remeshbottom = row number in grid (counting from the top),
                     above which the vertical grid spacing is
                     dilated or compressed to accomodate the
                     perturbations specified; at this row and below,
                     the grid spacing is not affected.
                     
                     NOTE: make sure you choose remeshbottom 
                     to be further down than the lowest point
                     of your defined surface topography.


CRPIT section
-------------

Defines the CRPIT parameters.  CRPIT=Controlled Remeshing with Partial
Interface Tracking.  See CRPIT notes for a more detailed explanation.
The CRPIT notes explain the concept of controlled remeshing and partial
interface tracking, the meaning of the control function and the definitions
of Regions 1, 2, 3, and 4.

This section starts with the word CRPIT on a line by itself.
The next lines are as follows:

        read(plus_file,*) crpit_interface, stablax_interface
           crpit_interface = number of the tracked interface to which
                             you want to apply CRPIT
           stablax_interface = lax stabilization parameter to use
                               for the CRPIT interface.  If you 
                               don't want to use any lax stabilization,
                               use stablax_interface = 0.d0
        
        read(plus_file,*) crpit_x1, crpit_x2, crpit_x3
           crpit_x1 = x value which defines the boundary between
                      Region 1 and Region 2 of the CRPIT control function
                      at time=0.
           crpit_x2 = x value which defines the boundary between 
                      Region 2 and Region 3 of the CRPIT control
                      function at time=0.
           crpit_x3 = x value which defines the boundary between 
                      Region 3 and Region 4 of the CRPIT control
                      function at time=0.               
        
        read(plus_file,*) crpit_f_fraction
           crpit_f_fraction = the fractional value of the CRPIT control
                              function in Region 1.
        
        read(plus_file,*) crpit_f_velocity
           crpit_f_velocity = the velocity at which the CRPIT control
                              function is moving
        
        read(plus_file,*) crpit_target
           crpit_target = the target value which the tracked interface 
                         (crpit_interface) moves toward, in Region 1.




File: surfaceproc_i

Defines input parameters for surface processes: erosion, sedimentation,
progradation, etc.

ERODE section
-------------

This section starts with the word ERODE (on a line by itself).
The next lines are as follows:

      read(76,*)erosion_model
           erosion_model = type of erosion model you want to use
                         = 0 for no erosion
                         = 1 for combinations of total erosion and
                             total deposition depending on
                             nerode_type
                         = 2 for slope-dependent erosion/deposition
                         = 3 for diffusive erosion (added by Susanne
                             Buiter)
                         
      read(76,*)erode_p(1)
        If erosion_model = 2 (slope dependent erosion) then
           erode_p(1) = erosion rate per second when the surface
                        slope = 1
        If erosion_model = 3 (diffusive erosion) then
           erode_p(1) = diffusion coefficient k (in m^2/s)
                  

      read(76,*)nerode_type, nerode_ends
         nerode_type = 1  erosion only
                     = 2  erosion and deposition (may be unstable, gives
                          element-by-element oscillations in Eulerian grid)
                     = 3  deposition only
         
         nerode_ends = 0  No erosion of surface endpoints
                          This is the recommended value when deformation
                          in the model will not propagate to the 
                          boundaries
                     = 1  Erodes points 1 and nx by same amount as
                          points 2 and nx-1, respectively.
                          This gives no-slope boundaries at the
                          edges for surface processes.
                     = 2  Erodes points 1 and nx using local slope
                          NOT CODED YET

      read(76,*)er_stablax, method_slope_avg
         er_stablax = Lax stabilization factor
                      The larger this value is, the more
                      stabilization is applied.
                        
         method_slope_avg = 1  calculate average slope by taking the 
                               absolute values of the slopes and then
                               averaging
                          = 2  calculate average slope by taking the
                               slope between points 1 and 3.
                               
         If this line (containing er_stablax and method_slope_avg) is
         missing from erode_i, the code should set er_stablax to 1
         and method_slope_avg to 1.  Look in the first few lines of
         the standard output (so.xout) to see what the code is using
         for er_stablax and method_slope_avg.
         

SEDIMENT section
----------------

This section starts with the word SEDIMENT (on a line by itself).

     read(76,*) nsedcolor
       nsedcolor = number of time/sediment material data lines
     
      do i=1,nsedcolor
        read(76,*) tchron_strat(i), sedcolor(i)
      enddo

      tchron_strat(i) = time (in seconds)
      sedcolor(i)     = the material number for the sediment that 
                        will be deposited when the model time is .ge. 
			tchron_strat(i).


PROGRADE section
----------------
This section starts with the word PROGRADE (on a line by itself).

      read(76,*) progradation_model

	  progradation_model = type of progradation to use 
			     = 0  no progradation
                             = 3  Lykke Gemmer's fill-and-spill
                                  (not implemented in code-feb23-04)
                             = 4  progradation + aggradation 

          NOTE: this version of SOPALE does not use the t_stop_prograde
          variable which was specified on this line in previous versions.
          Also, progradation models 1 and 2 are no longer recognized.
          Progradation model 4 is a replacement for progradation model 2.

            =======================================================
                       NOTES on Progradation Model 4

            This is a half-Gaussian progradation, plus uniform 
            aggradation, that vary with time according to a time-
            stepping function.

            The general form of the sediment profile function is

            h(x) = (h0(i)-hinf(i)) 
                      * exp(-((x-x0(i))**2)/prograde_length(i)**2))
                   + hinf(i) + haggrade(i) 

            The values of h0(i), hinf(i), x0(i), prograde_length(i), 
            and haggrade(i) are reset at the beginning of each time 
            interval(i) and do not vary during the time interval.

            During the time interval x0(t) and haggrade(t) vary 
            with time since the beginning of the time interval 
            according to

             x0(i, deltat) = x0(t=t(i)) + vprograde(i) * deltat
             haggrade(i,deltat) = haggrade(t=t(i)) + vaggrade(i)*deltat

            where 
               deltat is the time elapsed since the beginning of 
                      this time interval
               t(i)   is the absolute model time at the beginning 
                      of the time interval

            =======================================================

        IF progradation_model = 4 then

         read(76,*) numprograde

            numprograde = number of times at which the progradation and
                        aggradation parameters are specified.

         do i = 1, numprograde
           read(76,*) tprograde(i), h0(i), hinf(i), x0(i), 
                      prograde_length(i),vprograde(i), haggrade(i),
                      vaggrade(i)

          tprograde(i) = model time (in seconds) at which the profile 
                         values are set to these values 
     
          h0(i) = surface height at x.le.x0(i) (metres)
          hinf(i) = limit of surface height at x.eq.infinity (metres)
          x0(i) = initial x position where the half-Gaussian part of the
                  profile starts (metres)
          prograde_length(i) = Gaussian width of the profile as defined
                               above (metres)
          vprograde(i) = progradation velocity (m/s)
                        velocity at which progradation function moves in
                        positive x direction
                        (0.317097d-9 = 0.01m/year)
                        Negative values give retrogradation of the
                        profile but do not result in erosion of sediments
                        above the profile.
          haggrade(i) = base level for aggradation (metres)
          vaggrade(i) = aggradation velocity (m/s) which is uniform
                        for all x velocity at which sediments aggrade
                        (0.317097d-9 = 0.01m/year)
                        Negative values give degradation of the
                        profile but do not result in erosion of sediments
                        above the profile.

          NOTE: 1) All heights and positions are with respect to the origin
                   (x=0, y=0)
                2) All parameter values are the values for the time interval i
                   and are reset at the beginning of time interval i+1.
                3) The form of the profile function is quite flexible but it
                   is necessary to the user to calculate the initial values
                   for each time interval carefully if steps/jumps in the
                   positions of the profile are to be avoided.
                   For example, x0(i+1) = x0(i) + vprograde * (t(i+1)-t(i))
                   correctly positions the profile at the beginning of
                   interval(i+1) to the position it had at the end of
                   interval(i).
                4) Parameter values specified for tprograde(numprograde)
                   apply for all times larger than this value.
                   



WATERLOAD section
-----------------
This section starts with the word WATERLOAD.

waterloading  = 1 use nodal weights for water loading
              = 0  don't use water loading

refh_water, rho_water, stab_waterweight
    refh_water  = sea level (for code-sep17-03, assumed to be parallel 
		  with the x axis)
    rho_water   = density of water
    stab_waterweight = stabilization factor for the water weight
                       (versions code-dec04-03 and later)

NOTE: code-sep17-03 version has modified water loading to act as
      a pressure normal to the (land) surface.  This only works
      if the model is constructed in the x-z coordinate system
      such that the surface of the water (refh_water) is horizontal
      (i.e. parallel to the x coordinate axis) and that gravity acts
      in the minus z direction. 
      [ Will only work if ang_gravity is -90.]


File: matlab_i  (optional)

            write(*,*) 'Found file matlab_i; reading...'
            read(49,*)
            read(49,*)nout       ! interval for output timesteps
            read(49,*)
            read(49,*)icolor1    ! Eul materials output
            read(49,*)
            read(49,*)icolor2    ! Lag materials output
            read(49,*)
            read(49,*)iecord     ! Eul coordinates output
            read(49,*)
            read(49,*)ilcord     ! Lagrangian coordinates output
            read(49,*)
            read(49,*)istrain1   ! Eulerian strain output
            read(49,*)
            read(49,*)istrain2   ! Lagrangian strain output
            read(49,*)
            read(49,*)istrainr   ! Eulerian strain rate output
            read(49,*)
            read(49,*)itemper    ! Eul temperature output
            read(49,*)
            read(49,*)ievel      ! Eul velocities output
            read(49,*)
            read(49,*)ilvel      ! Lag velocities output
            read(49,*)
            read(49,*)isxx
            read(49,*)
            read(49,*)isyy
            read(49,*)
            read(49,*)isxy
            read(49,*)
            read(49,*)inodpres    ! Eul nodal pressures
            close(49)
            write(*,*) 'Closed file matlab_i'

Sopale documentation