Sopale documentation
USER GUIDE FOR SOPALE1 (Serial OPALE)
Version: code-nov26-03
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_pre_prs is specified as a material property
and should be chose 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 FLUI 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.
===========================================================================
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, t_stop_prograde
progradation_model = type of progradation to use
(0 = no progradation)
t_stop_prograde = time at which to stop progradation
if (progradation_model.eq.1.or.progradation_model.eq.2) then
read(76,*,IOSTAT=iostatus) h0,hinf,x0,prograde_length,vprograde,
& h0star ! [mod BL 2002/11/27]
h0 = surface height at x.eq.0 (metres)
hinf = limit of surface height at x.eq.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
All heights and positions are with respect to the origin
(x=0, y=0)
vprograde = progradation velocity (m/s)
velocity at which progradation function moves in
positive x direction
(0.317097d-9 = 0.01m/year)
h0star = lowered height behind the prograding surface
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
refh_water = sea level (for code-sep17-03, assumed to be parallel
with the x axis)
rho_water = density of water
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'