SOPALE
Sopale Nested
Processing SOPALE Output
Printing
MOZART
TMM (Thermo-mechanical Model)
old SOPALE documentation
Geodynamics home page
|
sopale_nested input
Notation
In many of the lines below, there are items in square brackets [].
The square brackets are meant to indicate that the item(s) inside
are optional. This DOES NOT APPLY to the section header line
[SOPALE1_i]. It must be entered as written here.
Lines whose first character is '!' are treated as comment lines;
they are flushed from the input stream.
[SOPALE1_i]
Indicates start of section.
template_i
template_i
number of lines in the model description, which follows immediately.
comment line
comment line
There should be template_i of these.
sopale_nested will echo these lines to the
standard output file. The modeller is encouraged to use these
lines as a description of the model.
nmodels
nmodels
sopale_nested reads and ignores this.
modelsn(1) modelsn(2) ... modelsn(nmodels)
modelsn(i), i=1,nmodels
sopale_nested reads and ignores this.
modelkey
modelkey
sopale_nested reads this with an a300 format, echos it to the
text output file, then forgets it.
amodel benchmark
amodel
sopale_nested reads and ignores amodel.
sopale_nested reads and ignores benchmark.
wres rres align_time [ update_viscosities [ cond_weight ]]
wres
= number of timesteps between writing of 2 restart files
rres
determines whether to read restart files.
= 0 means don't read restart files.
= 1 means read restart files.
align_time
(only has effect when rres = 1)
= 0 On restart, set the time to zero and
timestep to 1. This is appropriate when you are using restart
files which just establishes the state you want.
= 1 On restart, set the time and timestep to
those found in the restart files. This is good for restarting an
interrupted run.
Since version 2009-07-26_22.9: When restarting
with the pressure pump on (type_bc=300),
sopale_nested will stop if align_time is not 1.
Since version 2009-07-26_26.0: When restarting
with the pressure pump on (type_bc=300),
and pump_time greater than dt, sopale_nested will stop if
align_time is not 1.
Said another way:
- When the pressure pump is on
- The requirement use aligntime relates to the value of
pump-time. If pump_time is greater than the timestep length
aligntime must be used....Rajesh-type models. If pump_time is
less/equal to the timestep length or not defined the user has
the option to use aligntime or not.
- If the pressure-pump is off the user has the option to use
aligntime or not
Since version 2009-07-26_22.
update_viscosities
An optional parameter. Rationale
viscosities are normally updated on the second, and subsequent,
iterations of every time step.
= 0 means do not update viscosities.
= 1 means update viscosities. This is the default.
Since version 2009-07-26_27.0
cond_weight
An optional parameter.
Expected value is between 0.0 and 1.0; default is 0.5
Previous versions of sopale_nested implement
cond_weight=1.0
When a finite element calculation is done to compute the
initial temperature field (type_bct0=3),
and thermal conductivities vary
with temperature, the calculation sometimes does not converge,
and the error flips back and forth between two solutions, one
positive, the other negative . In an attempt to prevent this
"flipping", we use a weighted average of the thermal
conductivities from the current and previous iterations. The
formula is
cond_weight*current +
(1-cond_weight)*previous
There are test results
of a (previously) non-converging temperature field, where the
convergence error alternated between 21.xxx and -21.yyy degrees.
While convergence for cond_weight=0.5 in some
tests is very quick, there is a good possibility that this is
chance, as opposed to good design.
Convergence or non-convergence of the finite
element calculation will depend both on the initial temperature(s)
and the boundary conditions. Your mileage will vary.
nsave1
nsave1
number of saves for Eulerian grid
(Maximum number of saves is 200)
isave_1(1) isave_1(2) ... isave_1(nsave1)
(isave_1(i), i = 1, nsave1)
[only needed if nsave1 > 0 ]
timesteps at which to save the Eulerian grid
nsave2
nsave2
number of saves for Lagrangian grid
(Maximum number of saves is 200)
isave_2(1) isave_2(2) ... isave_2(nsave2)
(isave_2(i), i =1, nsave2)
[only needed if nsave2 > 0 ]
timesteps at which to save the Lagrangian grid
nameout1 [ nameout1_r ]
nameout1
soaple_nested uses this as a run name. Output file names are
prefixed with it. Copies of the input file(s) are made; the copies
have this pre-pended to their original name. Restart file names
may begin with it. The readability of file lists is improved if
the last character is an underscore.
length limit is 200.
nameout1_r
An optional input. run name for the restart files. '_i' will
be appended to this to create the input file names.
Default: nameout1
Since: implementation date not immediately at hand
time_nstep dt outunt [ dt_adjust dt_min_f dt_max_f dt_incr_f time_end [ dt_ref ]]
time_nstep
number of timesteps in this model run,
unless dt_adjust is 'yes'
If dt_adjust is 'yes' the model timestep lengths will vary
and the model will be run under cfl (Courant_Friedrichs_Lewy factor)
control to limit the maximum advection per timestep.
dt
timestep length in seconds, or initial timestep length if
dt_adjust is 'yes'
Restarts with dt_adjust = no
- dt (from this file) will be used.
Restarts with dt_adjust = yes
- If dt is found in the restart file, it will be used,
provided it fits within the range defined by dt, dt_min_f,
dt_max_f.
- If dt is found in the restart file, and does not fit within
the range defined by dt, dt_min_f, dt_max_f, then dt (from this
file) will be used.
outunt
sopale_nested read and ignores this.
dt_adjust
yes => The time step length (dt) will be adjusted each time
step
any other value => dt will not be adjusted.
First implementation version 2013-05-30_1.0 . frame numbers
on restart may be inconsistent.
retrofitted to version 2013-03-21_1.1 . frame numbers
on restart may be inconsistent.
Revised in version 2013-06-28_1.0 . frame numbers
on restart may be inconsistent.
Revised in version 2013-09-17_1.0
In version 2013-09-17_1.0:
dt_min_f
dt will not be allowed to become smaller than dt * dt_min_f.
only used if dt_adjust is 'yes'.
Default: 0.1
Suggested: 0.1
dt_max_f
dt will not be allowed to become larger than dt * dt_max_f.
only used if dt_adjust is 'yes'.
Default: 10.0
Suggested: 10.0
dt_incr_f
dt_incr_f moderates the rate at which dt is increased.
only used if dt_adjust is 'yes'.
Default: 0.1
Suggested: 0.1
time_end
When dt is adjustable, using the number of time steps to
determine the length of a model run is not useful. So this number
(in million years) is used.
Default: no default
dt_ref
When dt_adjust='yes', dt_ref is used to map timestep numbers
to time when writing output files. When doing restarts, setting
this number to the dt used to calculate the time step numbers for
each type of output file will ensure that frame numbers on the
output files behave consistently.
Default: dt
added in Sept 2013
bigu bigt
bigu
multiplier for penalizing velocity boundary conditions
Recommended value: 1d7
bigt
multiplier for penalizing the temperature boundary conditions
Recommended value: 1d7
length height nx1 ny1 betaxg betayg length2
height2 nx2 ny2
The variables in this line are a partial definition of the
eulerian and lagrangian grids.
length
length of Eulerian rectangle
height
height of Eulerian rectangle
nx1
number of columns of Eulerian grid
ny1
number of rows of Eulerian grid
betaxg
sopale_nested reads and ignores this.
betayg
sopale_nested reads and ignores this.
length2
length of Lagrangian rectangle
height2
height of Lagrangian rectangle
nx2
number of columns of Lagrangian grid
ny2
number of rows of Lagrangian grid
NOTE: surface is at "height", not 0
nodes are numbered starting at top surface
x=0 is at left edge
npert
npert
sopale_nested reads and ignores this.
nfpert
nfpert
sopale_nested reads and ignores this.
iformul
iformul
a parameter determining that plane strain calculation will used in Sopale
Required value = 2
hists(1) hists(2) hists(3) ... hists(11)
(hists(i),i=1,11)
sopale_nested reads and ignores this.
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, in SI units (m/s2
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
NOTE: [Jun 12, 2006] In versions of the code other than the ones listed below, there are two bugs
which prevent psmooth from having any effect (i.e. no pressure smoothing is done). The following versions of the code
have both these bugs fixed:
- code-feb07-06-adiabat, where the executable is dated later than Jun 8, 2006
- code-mar01-06-pttpaths, where the executable is dated later than Jun 8, 2006
See the email here for information on these bugs if you need to fix it in your
own code. The first email in the link deals with the first bug; the second email deals with the second bug, which was
discovered and fixed later.
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
see NOTE on the use of edotmax
n_matsets
n_matsets
number of material sets (aka rock types) whose mechanical
properties are to be defined
A material set is a collection of materials whose
properties are all the same. This allows the user to define
several materials which have the same mechanical properties
but different material numbers. This can be useful in
tracking particular areas during model evolution.
Said another way, a material set is a rock type with one
set of mechanical properties, but with several mineralologies.
The different mineralologies are represented by different
"colors".
NOTE: FOR VERSIONS OF THE CODE BEFORE code-nov16-05-matsets,
material sets were not used. Instead, the variable
color_number was used to define the
total number of materials.
mixhyp refsrate redphi nredphi plas_fluid_method
mixhyp
mixing model for viscous material properties (viscous
rheologies)
= 2, meaning viscous rheologies act in parallel.
Each viscous rheology is assumed to have the same second
invariant of 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.
Since version 2009-07-26_8:
= 1, meaning viscous rheologies act in series.
Each viscous rheology is assumed to have the same second
invariant of stress (they act in series). The total second
invariant of strain rate is the same as the strain rates of each
viscous rheology acting alone. The total viscous stress is
compared with the plastic yield stress to determine whether flow
is viscous or plastic.
This variable used to be floating point. Values like '2.0d0'
have been seen in older input files. The variable is now an
integer, and values like "2.0d0" should be changed to "2".
refsrate
reference strain rate to start viscosity computation viscous
creep flow models.
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)
NOTE: redphi and nredphi are part of a 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. See note on using redphi and
nredphi.
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 |
Description |
Comments |
| 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*sin(phi) + coh
|
|
1
Note: coh is not multiplied by cos(phi) |
Assume lambda_pore_prs = fluid pressure/solid pressure (mean stress).
If lambda_pore_prs <= 10**-30, then
splas = epressure*sin(phi) + coh
else
splas = (epressure - pwater)*(1-lambda_pore_prs)*sin(phi)+ coh
|
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.
Models using plas_fluid_method = 3 and plas_fluid_method =
4 may have problems converging in the first timestep
unless lambda_blend1 and lambda_blend2 are used. The
problem may be more pronounced with values of
lambda_pore_prs greater than 0.4. |
| 2 |
Assume lambda_pore_prs = fluid pressure/solid pressure (mean stress).
If lambda_pore_prs <= 10**-30, then
splas = epressure*sin(phi) + coh*cos(phi)
else
splas = (epressure - pwater)*(1-lambda_pore_prs)*sin(phi) + coh*cos(phi)
|
3
Note: coh is not multiplied by cos(phi) |
Assume
lambda_pore_prs = fluid pressure/overburden weight. In the code,
overburden weight is called plitho.
If lambda_pore_prs <= 10**-30, then
splas = epressure*sin(phi) + coh
else
splas = (epressure - (pwater + lambda_pore_prs*plitho))*sin(phi) + coh
|
| 4 |
Assume lambda_pore_prs = fluid pressure/overburden weight. In the code,
overburden weight is called plitho.
If lambda_pore_prs <= 10**-30, then
splas = epressure*sin(phi) + coh*cos(phi)
else
splas = (epressure - (pwater + lambda_pore_prs*plitho))*sin(phi) + coh*cos(phi)
|
| 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 - (pwater + lambda_pore_prs*phydro))*sin(phi) + coh*cos(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!
|
| Definitions |
splas is the magnitude
of the second invariant of the stress required for plastic
yielding.
epressure is the mean
stress calculated by the finite element solution. It
includes the effect of water loading because waterweight
is added to the body force vector in subroutine
"compute_body_forces".
pwateris the weight of the
water (rho_water*d*-gravity) applied as a vertical load in
"compute_body_forces".
plithois the weight
(lithostatic pressure) of the model computed at the centre
of the Eulerian elements. It DOES NOT include the water
load.
phydrois the weight of the
model assumed to be composed of water (rho=1000kg/m**3)
computed at the centre of the Eulerian elements. It DOES
NOT include the water load.
|
|
Material Sets
For each material set (i =
1, n_matsets) include these
lines of input.
color_numbers_for_material_set
color_numbers_for_material_set
This is a list of material numbers, indicating the material
numbers which belong to this material set. The items in the
list are separated by commas. Each item is a material number
or a range of material numbers. Example: 1,5,17-21
This indicates that materials 1, 5, 17,18,19,20,21 belong to
this material set.
density [ maximum_yield_stress [ mixing_reaction_scaling ]]
density
density of materials in material set i
maximum_yield_stress
maximum yield stress
Default: 1.0d12
Since: June 2014
mixing_reaction_scaling
mixing_reaction_scaling =1 optional flag only required if
this material set is subject to weakening that results from mixing
or reaction (metasomatism) in an arc_region. In such cases the
material density changes as a consequence of mixing or reaction
(see ARC section of User Guide). In this case the variable density
(above) changes as the model evolves and is interpreted as the
current density.
If (mixing_reaction_scaling. eq. 1) the evolving density is
used to vary the scaling factor nlcreep_sc in the nonlinear
(power-law) creep flow law (see input below) in which case the
scaling factor varies linearly from nlcreep_sc1 to nlcreep_sc2
over the range that the density changes. ( No scaling will occur
if nlcreep_sc1 =nlcreep_sc2.). The density scaling replaces the
strain-weakening scaling of nlcreep_sc.
The scaling for strain weakening is
nlcreep_sc = nlcreep_sc1 -
((straine – nlcreep_strain1) / (nlcreep-strain2 –nlcreep_strain1))
x (nlcreep_sc1 – nlcreep_sc2) Equation 1
In the case of density weakening ( mixing_reaction_scaling)
the scaling becomes
f_current = f_initial - ((rho3_current –
rho3_initial) / (rho3_final – rho3_initial)) x (f_initial –
f_final), Equation 2
so that the input parameters listed below are
reinterpreted. Note the current_strain (straine) is replaced by
the current_denisity (density)
If the density of this material set varies as a consequence of
mixing or reaction its density will be updated automatically by
the ARC module (if this is a restart the density = current_density
will be read from the restart file)
Default: 0
Since: May 2015
bulkv visc_on(1) visc_on(2) visc_on(3)
visc_on(4) [ bulkv_factor ]
bulkv
"bulk viscosity"
This is not actually a realistic value of bulk viscosity, but
rather a numerical parameter used to maintain incompressibility
of the system during solution by the penalty method. This value
is typically about 7 orders of magnitude larger than the lowest
viscosity in the model.
"visc_on(i), i=1,4" determines which of the viscous flow laws
are active during the calculation.
visc_on(1)
Flag for flow law number 1; grain size sensitive creep.
0 => do not use.
1 => use.
visc_on(2)
Flag for flow law number 2; nonlinear creep.
0 => do not use.
1 => use.
visc_on(3)
Flag for flow law number 3; Peierl's law creep.
0 => do not use.
1 => use.
visc_on(4)
Flag for flow law number 4; pressure solution creep.
0 => do not use.
1 => use.
bulkv_factor
bulkv_factor is a number less than or equal to 1.0d0 (the
default) used to multiply and rescale the bulk viscosity (bulkv) in
Eulerian elements with this mechanical material when there is an
enforced volume change as part of the mechanical finite element
solution (ie delv_vdt_phase and/or delv_vdt_arc are not equal to
0.0d0). The scaling scaling makes the element more compressible and
also reduces the nodal forces required to enforce the volume
change. It does not apply to delv_vdt_compact.
Default: 1.d0
Since: May 2015
phi1 phi2 phie1 phie2
lambda_pore_prs1 lambda_pore_prs2
lambda_pore_prse1 lambda_pore_prse2
phi1
angle of friction in degrees corresponding to strain level phie1
phi2
angle of friction in degrees corresponding to strain level phie2
phie1
strain level 1
phie2
strain level 2, which must be greater than phie1
Diagram:
________________ phi2
/|
/ |
/ |
/ |
phi1 ______________/ phie2 (strain level 2)
|
phie1 (strain level 1)
lambda_pore_prs1
pore pressure ratio corresponding to strain level lambda_pore_prse1
lambda_pore_prs2
pore pressure ratio corresponding to strain level lambda_pore_prse2
lambda_pore_prse1
first strain level
lambda_pore_prse2
second strain level, which must be greater than lambda_pore_prse1
________________ lambda_pore_prs2
/|
/ |
/ |
/ |
lambda_pore_prs1 __________/ lambda_pore_prse2
|
lambda_pore_prse1
More info: see note on use of phi and lambda_pore_prs
See note on plastic yield criterion and pore fluid pressure
coh1 coh2 cohe1 cohe2
coh1
cohesion corresponding to strain level cohe1
coh2
cohesion corresponding to strain level cohe2
cohe1
first strain level
cohe2
second strain level, which must be greater than cohe1
Diagram:
________________ coh2
/|
/ |
/ |
/ |
coh1 ______________/ cohe2
|
cohe1
GS gscreep_mu gscreep_m gscreep_T
gscreep_strain1 gscreep_strain2 gscreep_a0 gscreep_a1 gscreep_a2
gscreep_updown_T gscreep_VdivR
[GRAIN SIZE SENSITIVE CREEP]
gscreep_mu
viscosity at level of grainsize a0
To suppress gs_creep in parallel, take gs_creep_mu = 0
NOTE: See Philippe's note on viscosities
gscreep_m
[2] power of grainsize dependency
[m* in formula below]
gscreep_T
[3] temperature at which mu is (1/e) x that viscosity
at a0 and T
[T* in formula below]
gscreep_strain1
[4] S1* (see fig below)
gscreep_strain2
[5] S2* (see fig below)
NOTE: Always use S1* < S2*
gscreep_a0
[6] parameter a0* [ a0* in formula below ]
gscreep_a1
[7] parameter a1* [see fig below ]
gscreep_a2
[8] parameter a2* [see fig below ]
gscreep_updown_T
For case exp (- Td /T*), gscreep_updown_T = +1
For case exp (- T*/Td ) gscreep_updown_T = -1
[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]
gscreep_VdivR
Parameter for adding effect of activation volume V to the
creep laws for computing viscosities.
VdivR = V*/R, where V* is the activation volume and R
is the universal gas constant 8.3144 J/(K x mol).
Use 0.d0 if you don't want to have effect of activation volume.
PL nlcreep_eps nlcreep_sc1 nlcreep_sc2 nlcreep_n nlcreep_T nlcreep_strain1 nlcreep_strain2 nlcreep_updown_T nlcreep_VdivR [ peierls_threshold peierls_factor ]
[NONLINEAR, 'POWER-LAW', CREEP]
nlcreep_eps
[1] in the formula below = Aps = ½ x 3**((n+1)/2) x Auni
Where Aps (A invariant form) is the
experimental value of A, Auni, converted to the plane strain or
tensor invariant form, assuming that the experiment that
determined A was uniaxial deformation. If experiment was not
uniaxial the conversion formula must be recalculated
In the nonlinear (power-law) creep formulation used here
the second invariant of the stress (sigma2) is scaled by
nlcreep_sc. This scaling factor represents how the stress and
viscosity are modified by comparison with the standard flow law
(nlcreep_sc = constant) by the following factors
1) no scaling of the flow law nlcreep_sc = constant
(nlcreep_sc1 = nlcreep_sc2 = constant);
2) strain weakening in which nlcreep_sc varies linearly
from nlcreep_sc1 to nlcreep_sc2 over the range of second invariant
of the strain nlcreep_strain1 to nlcreep_strain2;
1) or 2) are assumed unless the input has the optional
input mixing_reaction_scaling = 1 (see density input)
3) weakening in which the above variables are interpreted
in terms of the effects of mixing or reaction in an arc region
controlled by mixing or reaction of this mechanical material
(mat3) with an added (injected) second material (mat2) that
progressively changes the density of this material as more of this
material is added (see description of mixing and reaction in ARC
section).
If (mixing_reaction_scaling. eq. 1) then this material set
is weakened by the effect of mixing or reaction such that
nlcreep_sc varies linearly from nlcreep_sc1 to nlcreep_sc2 over a
range of density defined by the inputs nlcreep_strain1 = density
(the density at which weakening begins) and nlcreep_strain2 is the
density beyond which no further density scaling occurs.
Note the inputs for the mechanical material properties
(listed here) may, or may not be consistent with the densities and
value of rho3_increase in the reaction part of the ARC
section. Consistency requires that nlcreep_strain2 =
initial_density and nlcreep_strain1 = initial_density +
rho3_increase. However, the user can define a range of density
increase that is different from the range defined in the ARC
section.
nlcreep_sc1
[2] in the formula below = initial strain weakening or if (
mixing_reaction_scaling . eq.1) initial scaling factor (
f_initial)
nlcreep_sc2
[3] in the formula below = final strain weakening or if (
mixing_reaction_scaling . eq.1) final scaling factor (
f_final)
(Note these can also be used to give strain, or
mixing_reaction_scaling hardening)
nlcreep_n
[4] in the formula below = power-law exponent n (ie value of
the power in the power law flow equation)
nlcreep_T
[5] in the formula below = T* = Q/R where Q is the activation
energy (kJ/mol) and R the universal gas constant, 8.3144 J/(K x
mol)
nlcreep_strain1
[6] in the formula below = value of strain for the onset of
strain weakening, or if ( mixing_reaction_scaling . eq.1) value
of density at which scaling starts (rho3_initial)
nlcreep_strain2
[7] in the formula below = value of strain at which strain
weakening ends, or if ( mixing_reaction_scaling . eq.1) value of
density at which scaling ends (rho3_final)
nlcreep_updown_T
[8] in the formula below = -1, 1, or 0 (see below) Use -1 for
normal calculations
nlcreep_VdivR
[for versions later than code-jul26-04-testV]
Parameter for adding effect of activation volume V to the
creep laws for computing viscosities.
[9] in equation above
VdivR = V*/R, where V* is the activation volume (SI units) and
R is the universal gas constant 8.3144 J/(K x mol).
Use 0.d0 if you don't want to have the effect of activation
volume.
peierls_threshold
If the calculated sigma2 is larger than this, it is multiplied
by peierls_factor. After the multiplication, sigma2
is forced to be at least peierls_threshold.
Since: July 2014
Restrictions: must be positive.
Default: 1.d20
peierls_factor
Since: July 2014
Restrictions: must be positive, and less than or equal to 1.d0
Default: 1.d0
The formula below is the flow law written as strain rate (edot)
calculated from stress (sigma) where both are the square root of the
second invariants.
[1] [4] [9]
| n* |
edot = Aps x (sigma/sc) exp –((T* + P x VdivR)/ T)
| |
sc1=[2] [5] if updown= -1
sc2=[3]
[2] sc1* _____________
|\
| \
| \
| \
[3] sc2* _ _ _ _ _ _ | _ _\__________________
| |
| |
S1* S2*
[6] [7]
P is the lithostatic pressure (Pa) and T is temperature (K) both
calculated by Sopale.
The actual formulae used in Sopale are the inverse of the
equation above, ie the stress Sigma2 is calculated from the strain
rate (square root second invariant). See subroutine viscosities1_2
if (nlcreep_updown_T.eq.-1) then
sigma2 = nlcreep_sc * ((edot/nlcreep_eps)**(1.d0/nlcreep_n)) *
dexp((nlcreep_T + plitho * nlcreep_vdivr) / (nlcreep_n * tlg(l)))
endif
if (nlcreep_updown_T.eq.1) then
sigma2 = nlcreep_sc * ((edot/nlcreep_eps)**(1.d0/nlcreep_n)) *
dexp((tlg(l))/(nlcreep_n * (nlcreep_T + plitho * nlcreep_vdivr))
endif
if (nlcreep_updown_T.eq.0) then
sigma2 = nlcreep_sc * ((edot/nlcreep_eps)**(1.d0/nlcreep_n))
endif
PL pei_eps pei_Tn pei_sc pei_T2 pei_strain1 pei_strain2
[PEIERL'S LAW CREEP]
pei_eps
[1] in diagram below
reference strain rate
pei_Tn
[2] in diagram below
[use 0 to suppress in parallel]
pei_sc
[3] in diagram below
pei_T2
[4] in diagram below
pei_strain1
not used
pei_strain2
not used
PEIERL'S LAW POWER-LAW BREAKDOWN
NOTE: E = epsilon
[2]
. . Tn*/T [4]
E = E* (sigma/sc*) exp ( - T2* / T )
[1] [3]
PS eps p d sc1 sc2 t updown_t strain1 strain2
PRESSURE SOLUTION CREEP
eps
is -
p
is -
d
is -
sc1
is -
sc2
is -
t
is -
updown_t
is -
strain1
is -
strain2
is -
. .
( e / eps ) = (sigma/sc) * exp -(T*/ T ) * (1/ (T . dp ) )
|
sc1 |____________
| |\
| | \
| | \
| | \
sc2 |_ _ _ _ _ _ | _ _\__________________
| | |
| | |
+------------+----+---------------------------
strain1 strain2
texpand tdensity melt_weaken temp_solidus
temp_melt_weakened visc_melt_weakened [ visc_bingham ]
texpand
volume coefficient of thermal expansion.
See notes on
Mantle Thermal Expansivity.
Since version 2013-02-08_1.0: If negative, the value
will be ignored and 4 parameters for a variable coefficient of
thermal expansion will be read from the next line.
temperature1 temperature2 alpha1 alpha2
temperature1
For temperatures <= temperature1, the coefficient
of thermal expansion will be alpha1.
temperature2
For temperatures between temperature1 and
temperature2, the coefficient of thermal expansion will
be a linear interpolation between alpha1 and
alpha2.
For temperatures => temperature2, the coefficient
of thermal expansion will be alpha2.
tdensity
reference temperature for density (in absolute degrees)
density is
rho(T) = rho0 (1 - alpha(T - tdensity))
melt_weaken
'0' => no melt weakening
'1' => material melt weakens
'2' => material melt weakens with extra parameters
Since June 2014: if melt_weaken=2, the extra parameters
will be read from the next line.
The pressure and temperature pairs (0, melt2_t0), (melt2_p1,
melt2_t1), (melt2_p2, melt2_t2) define a line in P/T space. That
line is the solidus of the material. The liquidus is solidus +
melt2_t_int.
melt2_t0 melt2_t1 melt2_p1 melt2_t2
melt2_p2 melt2_t_int
temp_solidus
temperature (degrees Kelvin) at which melt weakening starts
temp_melt_weakened
temperature (degrees Kelvin) at which melt weakening is complete
visc_melt_weakened
linear viscosity of melt-weakened material
visc_bingham
bingham viscosity. Since version 2009-07-26_22.2. Optional. Default is 0.0
END OF INPUT FOR EACH MATERIAL'S MECHANICAL PROPERTIES
type_bc
type_bc
Defines the type of boundary condition to apply to the model.
The following table summarizes the type_bc values. For
detailed descriptions of the type_bc input, please see the
document on "SOPALE
Boundary Condition Input"
type_bc |
Input variables and summary description |
300 |
vxleft_up vxleft_down yleft1 yleft2
rollleft vxright_up vxright_down yright1 yright2 rollright
vytop0 vxbottom0 reference_pressure_column
reference_material [ pump_time [ isos_time_scale
[ delta_vel_factor [ delta_vel_toler [ start_time flag ]]]]]]
The recommended, updated version of
type_bc = 3.
Defines horizontal velocities at side boundaries (can have
inflow and outflow); can set top boundary to rigid top (no
vertical movement); can set bottom boundary to no-slip (no
horizontal movement)
Other features: side roller boundaries are more
consistently implemented than in type_bc=3; "pressure pump"
feature adjusts horizontal velocities at side boundaries to
keep the pressure within the model constant. Details for
type_bc = 300.
Since version 2009-07-26_21.2 an optional parameter:
pump_time
Since version 2009-07-26_23.0 an optional parameter:
isos_time_scale
Since version 2009-07-26_29.0 an optional parameter:
delta_vel_factor
Since version 2012-06-08_2.0 an optional parameter:
delta_vel_toler. See type_bc 301 for details.
Since July 2014; two optional parameters
start_time
If present, sets the start time for the boundary condition
definition on this line.
flag
'more' => indicates there will be another boundary
condition. Spelling flexible provided the word is not 'last'.
'last' => indicates this is the last boundary condition.
|
|
301
(current)
|
(Since version 2013-01-29_1.0)
An extension of type 300.
In type 300, constant velocities are applied to the top and
bottom of each side of the model. The velocities are choosen
so that the mass (or volume, if density is assumed to be
constant) in the model remains constant.
In type 301, the velocities applied to the top of each
side of the model are adjusted every time step so that the
forces applied approach a 'target_force'. The bottom
velocities are then adjusted so that the mass (or volume, if
density is assumed to be constant) in the model remains
constant. Forces are measured from the top down to
yleft1 on the left side, and from the top down
to yright1 on the right side.
Details
for type_bc=301
vxleft_up target_force_left yleft1 yleft2 rollleft
vxleft_up
Initial value of vxleft_up (like bc_type 300)
target_force_left
target force (from the top down to yleft1) for the left side.
'ft' in the calculation
of vxleft_up
(others like bc_type 300)
vxright_up target_force_right yright1 yright2 rollright
vxright_up
Initial value of vxright_up (like bc_type 300)
target_force_right
target force (from the top down to yright1) for the right side.
'ft' in the calculation
of vxright_up
(others like bc_type 300)
vytop0 vxbottom0
(like bc_type 300)
reference_pressure_column
reference_material [ pump_time [ isos_time_scale [ delta_vel_factor
[ delta_vel_toler ]]]]
The parameters in square brackets are optional
after version 2013-01-29_1.0. In version
2013-01-29_1.0 they are mandatory.
reference_pressure_column
(like bc_type 300)
reference_material
(like bc_type 300)
pump_time
(like bc_type 300)
isos_time_scale
(like bc_type 300)
delta_vel_factor
(like bc_type 300)
delta_vel_toler
When the pressure pump calculates delta_vel, the
calculation includes subtracting plithob in the reference
column from the average plithob. If abs ((ref_plithob -
ave_plithob) - 1.0) < delta_vel_toler, the calculation is
skipped and delta_vel is set to zero.
mantle_wind_vel fraction_left
mantle_wind_vel
mantle_wind_vel (in the usual SI units). This velocity
is added to the velocites in the lower region ....below the
lithosphere, after all of the other velocities have been set
using the rescaling velocity. Note that the addition of this
velocity is with the same sign on both sides....in one side,
out the other.
fraction_left
Determines how the bottom velocities are distributed
between the left and right sides of the model.
'separate': Each side will be balanced separately.
After vxleft_up is re-calculated, vxleft_down will be set so
that volume_in = volume_out on the left side. The right
side calculation will be similar.
<number>: The volume change due to the upper
velocities will be distributed (in the lower region) between
the left and right sides of the model.
- fraction_left * volume_change will go into, or out of,
the left side.
- (1 - fraction_left) * volume_change will go into, or out
of, the right side.
0.5 will give equal volumes on each side.
use_tectonic_forces force_column force_process [ filter_length filter_alpha ]
use_tectonic_forces
= 1: use tectonic forces.
force_tect = total horizontal driving force (force_dyna) - mean stress (force_epress).
See .out file for values of force_tect,
force_dyna and force_epress
= 0: use total forces
force_column
column number (measured from both left and right side) to
be used when calculating the forces. Column 1 may give very
jittery results. Column 10 is suggested.
force_process
After the force type (use_tectonic_forces)
and the column number (force_column) have been
selected, there is a choice of how to process the result
before using it in the velocity rescaling formula.
'raw': use the force calculated from the current time step.
'filtered': use the output of a low-pass time series
filter, where the time series is the boundary force of the
previous filter_length time steps. The filter
parameter filter_alpha can be used to set the filter
sensitivity.
filter_length
Default: 20
filter_alpha
Initial tests suggest that a value 0.2 gives
some lag, and much smoothing, whereas 0.4 gives less lag and
less smoothing.
Default: 0.4
vxleft_up_min vxleft_up_max vxleft_up_velocity_scale
vxleft_up_min, vxleft_up_max
In each timestep after vxleft_up is
re-calculated, the result is compared to these minimum and
maximum velocities and set to one of these values if
'out of range'.
vxleft_up_velocity_scale {in version 2013-01-29_1.0: vxleft_up_factor}
'Vd' in the calculation of vxleft_up.
vxright_up_min vxright_up_max vxright_up_velocity_scale
vxright_up_min, vxright_up_max
In each timestep after vxright_up is
re-calculated, the result is are compared to these minimum and
maximum velocities and set to one of these values if
'out of range'.
vxright_up_velocity_scale {in version 2013-01-29_1.0: vxright_up_factor}
'Vd' in the calculation of vxright_up.
|
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.
On a restart, if align_time=0, smoothing
will be done if ismooth is > 0. 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 and later, 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
Recommended value = 1
leftsmooth
Eulerian column number defining the left edge of the portion of the surface to be smoothed.
The minimum value is 2.
rightsmooth
Eulerian column number defining the right edge of the portion of the surface to be smoothed.
The maximum value is nx - 1.
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 ≤ verror
where maxdv is the greatest change in velocity
(in either x or y direction) since the previous iteration
verror
relative error for timesteps greater than 1 in Picard iteration
verror1
relative error for timestep 1 in Picard iterations
miterv
maximum 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.
Actually, the code always does one more iteration
than miterv. This could be called a bug. However
the behaviour (bug) has existed for so many years that it seems the
lesser of two evils to leave the code as is.
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
maximum number of iterations per timestep before the viscosity
is reinitialized and computation is started from scratch.
Should always check by setting kinitvis = 1
The above is the original documentation. At present
(2010-04-16) kinitvis is a one-shot parameter; viscosities will be
initialized when the iteration number = kinitvis, in addition to
when they would otherwise be intialized.
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 beyond the first timestep
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 minumum number of iterations
(normal behaviour)
force_continue
This is a flag to indicate 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.
ncyc2 color2op lagtest eulcolorrule eul_density_rule
ncyc2
number of cycles in Lag/Eul remeshing/tracking per step
Suggested value = 2
color2op
sopale_nested reads and ignores this.
lagtest
sopale_nested reads and ignores this.
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 plurality, then the Eulerian cell's existing color
is not changed.
eul_density_rule
1 = calculate cell mass using the densities of its Lagrangian
particles.
Note that each Lagrangian particle has two
densities; one used in the finite element calculations, the other
calculated by the compaction module.
2 = calculate cell mass using eulerian cell material densities. Obsolete.
Must be 1 if there are any material phase change rules.
Since 2016-06-30: Must be 1.
ucompute tcompute ale_compute dt_compute tbottom
ttop e_temperature cfl lambda_blend1 lambda_blend2
ucompute
sopale_nested stops if this is not 1.
tcompute
1 = compute a temperature field using an FE computation
0 = no FE computation
Use 0 if you are running a mechanical-only model (no thermal-mechanical coupling)
ale_compute
sopale_nested stops if this is not 1.
dtcompute
sopale_nested reads and ignores this.
tbottom
temperature at bottom boundary (in Kelvin)
ttop
temperature at surface (in Kelvin)
e_temperature
2014-02-03: disabled. sopale_nested reads and ignores this.
cfl
Courant–Friedrichs–Lewy factor; used to adjust timestep length
(dt) when dt_adjust = 'yes'. suggest approx 0.2, no greater than 0.3
lambda_blend1
Iteration number during the first timestep where the
calculation of plastic yielding starts to incorporate the value
of lambda_pore_prs (fluid
pressure ratio)
lambda_blend2
Iteration number during the first timestep where the
calculation of plastic yielding has finished incorporating the
full value of lambda_pore_prs
(fluid pressure ratio).
lambda_blend1 and lambda_blend2 define an 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 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.
steady analy esteady
steady
sopale_nested reads and ignores this.
analy
sopale_nested reads and ignores this.
esteady
sopale_nested reads and ignores this.
OUTFIL
OUTFIL
sopale_nested reads and ignores this.
MATFIL
MATFIL
sopale_nested reads and ignores this.
ICASE
ICASE
sopale_nested reads and ignores this.
CACHSZ
CACHSZ
sopale_nested reads and ignores this.
LEVEL
LEVEL
sopale_nested reads and ignores this.
ibinread
ibinread
sopale_nested reads and ignores this.
filename1
filename1
sopale_nested reads and ignores this.
filename2
filename2
sopale_nested reads and ignores this.
filenameg1
filenameg1
sopale_nested reads and ignores this.
filenameg3
filenameg3
sopale_nested reads and ignores this.
filenameg4
filenameg4
sopale_nested reads and ignores this.
mechanical color assignment
Here we assign a initial mechanical color (not material) to the
lagrangian particles.
ncolor2 [ s1 s2 ... ]
ncolor2
number of Lagrangian mechanical color boxes defined on Lagrangian grid
s1 s2 ...
seed used to initialize the pseudo-random number generator. On
the scrXm family and the stable, when using binaries compiled with
Intel's ifort, this should be two integers, maximum 9 digits each.
Other compilers may need a different number. gfortran's
pseudo-random number uses a seed with eight integers. Eight is the
maximum number allowed for at present.
Only used if strain2_type is 'random'.
Default: 123456789 362436069 521288629 316191069 987654321 458629013 582859209 438195021
Since version 2009-07-26_29.0
For each box (i = 1, ncolor2), there are 3 lines of input.
Coordinates are in meters, time is in seconds.
x1 y1 x2 y2
x1 y1
coordinates of upper left corner
x2 y2
coordinates of lower left corner
x3 y3 x4 y4
x3 y3
coordinates of lower right corner
x4 y4
coordinates of upper right corner
newcol2 [ strain2_type strain2_magnitude [ runtime-flag time ]]
newcol2
the mechanical color number assigned to all Lagrangian particles in this box
strain2_type
'constant: The strain of the Lagrangian particles will
be initialized strain2_magnitude.
'random': The strain of the Lagrangian particles will be
initialized to a number from the pseudo-random number generator in
the range -strain2_magnitude to
strain2_magnitude.
Default: 'constant'
Since version 2009-07-26_29.0
strain2_magnitude
magnitude used when initializing the strain of Lagrangian particles.
Default: 0.d0
Since version 2009-07-26_29.0
runtime-flag
'F' => the box will be applied at model initialization.
'T' => the box will be applied at time='time'
Default: F
Since version 2014-03-??
time
Time (in seconds) at which the box will be applied. Only used if
'runtime-flag' is 'T'
Default: 0.0
Since version 2014-03-??
thermal material assignment
Here we assign a initial thermal material to the lagrangian
particles.
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).
ncolor2
ncolor2
number of thermal material boxes defined on Lagrangian grid
For each box (i = 1, ncolor2), there are 3 lines of input.
Coordinates are in meters, time is in seconds.
x1 y1 x2 y2
x1 y1
coordinates of upper left corner
x2 y2
coordinates of lower left corner
x3 y3 x4 y4
x3 y3
coordinates of lower right corner
x4 y4
coordinates of upper right corner
newcol2 [ runtime-flag time ]
newcol2
the thermal material number assigned to all Lagrangian
particles in this box
runtime-flag
'F' => the box will be applied at model initialization.
'T' => the box will be applied at time='time'
Default: F
Since version 2014-03-??
time
Time (in seconds) at which the box will be applied. Only used if
'runtime-flag' is 'T'
Default: 0.0
Since version 2014-03-??
dmtm n_skip_thermal
dmtm
should be 1
n_skip_thermal
number of timesteps 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 timesteps)
Only if dmtm = 1, then read the following:
type_bct0 [ temperature_initial maxter_i [ max_th_iters ] ]
type_bct0
Defines how the initial Eulerian temperature field is
calculated.
0 => the initial temperature field is computed based on the
inputs in the ETEMP
section of
sopale1_plus_i. This gives
a 2D initial steady-state temperature variation with
depth.
3 => use a finite element calculation to compute the initial
steady-state temperature field.
Note that
the ETEMP values for temperature will be overwritten if
type_bct0=3.
the type_bct0=3 temperatures will be overwritten if EPERTURB-T is present.
For more about temperature initialization, see here.
Since version 2009-07-26_10
temperature_initial
initial temperature for the steady state thermal solution.
Since version 2009-07-26_22: If this is negative, the
ETEMP field values are used as the initial temperature for the
steady state thermal calculation.
The default is -1.0
maxter_i
The steady state thermal solution is declared converged when
the maximum change in temperature is below this.
Since version 2009-07-26_22: Default is 1.d0 .
max_th_iters
maximum number of thermal iterations. Used to limit both the
initial steady state and time-stepping thermal solution. Default
value is 100.
Only If type_bct0=3, the next line
defines parameters for the boundary conditions to use
when solving for the initial steady-state temperature
field.
bct30(1) bct30(2) bct30(3) bct30(4) bct30(5)
bct30(6) bct30(7) bct30(8) bct30(9) bct30(10) bct30(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)
Since version 2009-07-26_10:
= 2 bottom is imposed temperature that varies with space. Read
extra input lines below.
= 3 bottom is imposed flux that varies with space. Read extra
input lines below.
bct30(3)
= temperature or flux at base (depending on value of bct30(2).
Since version 2009-07-26_10:
If bct30(2) = 2 or 3 the temperature or flux varies in space, and
bct30(3) is ignored. Read input below.
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
Only if bct30(2) = 2 or 3
Since version 2009-07-26_10
Spacial variation of basal temperature (bct30(2) = 2) or heat flux (bct30(2) = 3)
num_space_points
= number of points in the piecewise function defining the
spacial variation of the basal temperature or flux.
x_base_thermal_positions(i), i=1,num_space_points
= x positions of piecewise function defining spacial
variation of the basal temperature or flux.
val_thermal_bc(i), i=1,num_space_points
= values of temperature or flux in the piecewise function.
ntperts
ntperts
number of boxes for temperature perturbation
if zero, no boxes (below) will be read.
For each box, read in coordinates of corners and temperature
perturbation value (tpert). The value of tpert is added to the
initial temperature values.
For more about temperature initialization, see here.
xx1 xx2 xx3 xx4
xx1
x coordinate of top left corner of box (units: metres)
xx2
z coordinate of top left corner of box
xx3
x coordinate of bottom left corner
xx4
z coordinate of bottom left corner
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
tpert
tpert
Add this to the temperatures inside the box.
This page was last modified on Friday, 01-Jul-2016 13:50:40 ADT
Comments to geodynam at dal dot ca
|