INPUT FILES FOR THE THERMO-MECHANICAL CODE
Created: July 21, 1998
Last modified: July 25, 2002 by Bonny Lee
General notes:
Since the thermo-mechanical code was based on the microfem code, there
are input files which may have been used in microfem but aren't useful
or relevant in the thermo-mechanical model. However, I haven't yet
gotten around to cleaning up the thermo-mechanical code and figuring out
which files really aren't needed. I've noted the "possibly irrelevant"
files below.
==============
File: Export_i
Unit: 247
Purpose: Possibly irrelevant to the thermo-mechanical code. The microfem
documentation says it contains one value, iexport_what. If
iexport_what = 1, then microfem will open a file "Export_s"
and write the u0 array after convergence.
==============
File: Import_i
Unit: 218
Purpose: Possibly irrelevant to the thermo-mechanical code? The microfem
documentation says it contains one value, import_what. If
import_what=1, then read imported vectors from Import_s
and compare with solution.
==========
File: NL_i
Purpose: Contains various parameters
Format/Content:
All lines are free-format.
Line 1: nlflag nl_read
nlflag = ?
nl_read = ? (doesn't appear to be used)
Line 2: itmaxn er0n igsm
itmaxn = maximum number of iterations per timestep
er0n = max error value for convergence
igsm = parameter for smoothing grid?
doesn't seem to be used in code (use is commented out)
Line 3: pcohes imix dmatx2
pcohes = cohesion pressure (in pascals)
imix = mixing scheme; 1= averaging, 2= weakest material dominates,
3=strongest material dominates
dmatx2 = ?
Line 4: ichange_mat ichange_geom ichange_bc ichange_tem
These appear to be flags which you can use to tell the code
to read in new material properties (ichange_mat), new boundary
conditions (ichange_bc), or new temperatures (ichange_tem).
I don't know if this is currently properly implemented (I have
the feeling it's not..) Currently, they're usually all set to 0.
Example:
0 1
120 1. 1 !#iter err ? (overwrites c_i)
5.e+07 2 1.e+18 !pcohes imix dmatx2 - overwrites values in g0_s
0 0 0 0
=========
File: X_i
Purpose: This is a control file for running microfem models interactively
and looking at outputs along the way. If you're running a job on the sp2
then you don't want this, so set i_X to a very large number (> max # ts).
When running jobs on one of the workstations, you sometimes like to see
output as the job progresses. Then you can use this file to control which
timesteps get plotted. Or you can specify that the plots are to be saved
directly to a file (Paula's modifications).
Format/Content:
All lines free-format
Line 1: i_X isitX ipsplot iXdebug
i_X - specifies plotting interval: plots are displayed at
mod(kstep-1,i_X)
isitX - same as i_X but for iterations
ipsplot=0 for plot to screen
=1 for Postscript output only
iXdebug - if iXdebug=1, then do a few more plots for debugging
Line 2: psout
psout = name of file to save Postscript output to, if ipsplot=1
(only read if ipsplot=1)
Example:
9000 200 0 0
p.out
==============================
File: check_numerical_errors_i
Unit: 220
Purpose: ?
Format/Content:
All lines free-format.
Line 1: ichecknl
ichecknl =1 save information to check_numerical_errors_s
=0 don't save
Line 2: ssy1 ssy2 ssy3
Line 3: visc_s pres_s srat_s
Example:
1 (1=save information to check_numerical_errors_s, 0=don't save)
1.0 1.05 1.1
1.e+00 1.e+00 1.e+00
===============
File: initial_thermal_parms_i
Unit: 28
Purpose: Defines parameters required for thermal preprocessor, which
sets up initial thermal conditions. This file replaces the
interactive preprocessing section in the 1998 version of
the thermo-mechanical code.
Format/Content:
Line 1: xlamused ismootop
xlamused = bulk modulus in Pa.s
ismootop = timestep interval for surface smoothing
Line 2: iupsla
Line 3: cp rau conduc
cp = heat capacity
rau = density ("rau" is Philippe's version of "rho")
conduc = conductivity
Line 4: t_peclet rlength
t_peclet = Peclet number
rlength = length scale
Line 5: fead_bc(2),fead_bc(3),fead_bc(4)
These are thermal boundary conditions
fead_bc(2) = left heat flux in W/m**2
fead_bc(3) = right heat flux in W/m**2
fead_bc(4) = bottom heat flux in W/m**2
Line 6: fead_bc1(4)
fead_bc1(4) = bottom right corner initial temperature (in
degrees Celsius)
Line 7: t_bath
t_bath = bath temperature
Line 8: isingu
isingu = column number + 1 of singularity in mechanical grid
(first column where basal velocity > 0)
NOTE: Cross-check with input in microfem_g0_i
Ensure the convergence velocity is constant and decreases to 0
over 1 element - this constraint doesn't exist for purely
mechanical runs.
Line 9: i_th_colors
i_th_colors = the number of thermal material types, in therset_E_2_i
[Check therset_L_1_i also]
Line 10:(i_modulate(kk),kk=1,i_th_colors)
i_modulate(kk) = 0 means don't modulate the nominal source
strength by the value of the time
function at time 0 for this thermal
material
= 1 means DO modulate the nominal source
strength by the value of the time
function at time 0 for this thermal
material
Example:
1e+32 80 (bulk modulus in Pa.s; frequency in ts for surface smoothing)
56 (iupsla)
750. 3000. 2.00 (heat capacity; density; conductivity)
0 35000. (Peclet number; length scale)
0. 0. 0.02 (Thermal boundary conditions: left, right, bottom heat fluxes in W/m2))
1313.83 (bottom right corner initial temperature)
1350. (bath temperature)
122 (isingu)
3 (material thermal colors)
1 1 1 (modulation array for each of the initial source (radioactive) strenghs in therset_E_2_radio_i and therset_L_1_radio_i)
===============
File: mecprop_i
Unit: 348
Purpose: Defines mechanical properties for a set of materials.
Format/Content:
Line 1: mset
Format: i2
mset = number of materials being defined
Lines 2 to (mset*2)+1:
Format: 6e9.2
For each material, we have 3 lines in the input file:
dmat(1) dmat(2) dmat(3) dmat(4) dmat(5) dmat(6)
dmat(7) dmat(8) dmat(9) dmat(10) dmat(11) dmat(12)
dmat(13) dmat(14) dmat(15) dmat(16) dmat(17) dmat(18)
where
dmat(1) = mu_infinity = upper bound on viscosity
dmat(2) = K = bulk modulus
dmat(3) = phi = angle of friction, in degrees
dmat(4) = cohesion
dmat(5) = density
dmat(6) = maximum yield stress (this acts as a filter)
dmat(7) = B*
dmat(8) = activation energy (in Joules)
dmat(9) = power law exponent
dmat(10) = approximate guess at epsilon dot, used in initial calculation
dmat(11) = not used
dmat(12) = not used
dmat(13) = temp_s = solidus temperature in degrees C
dmat(14) = delta_s_l = melting range, temperature interval over which
viscosity decreases during melting
dmat(15) = visc_l = viscosity of material for temperatures equal
and above temp_s + delta_s_l
dmat(16) = temp_s_rau = temp at which density starts to decrease on
melting
dmat(17) = delta_s_l_rau = melting range, temperature interval over
which density starts to decrease on melting
dmat(18) = delrau = maximum decrease in density during melting
in kgm/m**3
The materials must be ordered starting with the strongest and decreasing
in strength.
NOTE: Versions of the code before July 22, 1998 had only 2 lines (i.e.12
properties) per material, of which dmat(11) and dmat(12) were not used.
Versions of the code after July 22, 1998 but before October 7, 1998
had only 2 lines (i.e. 12 properties) per material, of which dmat(11)
was visc_l and dmat(12) was delta_s_l.
Code version TMO.2.0 is the first to use 3 lines (18 properties) per
material.
Example:
6 mat priority high 10 low 1 therefore list strongest first, weakest last
1.00e+28 0.10e+19 1.50e+01 1.00e+00 3.30e+03 1.00e+38
7.75e+06 4.98e+05 4.48e+00 1.00e-19 0.00e+00 0.00e+00 This is mantle b*ol
1.00e+28 0.10e+19 1.50e+01 1.00e+00 2.70e+03 1.00e+38
1.05e+05 3.35e+05 2.60e+00 1.00e-19 0.00e+00 0.00e+00 This is lower crust b*pyx
1.00e+28 0.10e+19 1.50e+01 1.00e+00 2.70e+03 1.00e+38
1.05e+05 3.35e+05 2.60e+00 1.00e-19 0.00e+00 0.00e+00 This is lower crust b*pyx
1.00e+28 0.10e+19 1.50e+01 1.00e+00 2.70e+03 1.00e+38
2.98e+06 2.39e+05 3.20e+00 1.00e-19 0.00e+00 0.00e+00 This is feldspar b*fspar
1.00e+28 0.10e+19 1.50e+01 1.00e+00 2.70e+03 1.00e+38
1.28e+08 1.51e+05 1.80e+00 1.00e-19 0.00e+00 0.00e+00 This is strong quartz b*qz*10
1.28e+07 1.51e+05 1.80e+00 1.00e-19 0.00e+00 0.00e+00 This is quartz b*qz
1.00e+28 0.10e+19 1.50e+01 1.00e+00 2.70e+03 1.00e+38
number of mechanical set numbers
for each enter 12 numbers in 6e.. format
upper viscosity bulk modulus angle of friction ... same as regular previously in microfem_g0_i
for imix=2, list mec props in order of strongest first, weakest last
(preference given to weak materials)
for imix=3, list mec props in order of strongest first, weakest last
(preference given to strongest materials)
imix is defined in NL_i
==================
File: mecset_E_1_i
Unit: 344
Purpose: Defines the placement of the different materials (which were defined
in mecprop_i) on the Eulerian grid.
Format/Content:
Line 1: numpro blowl
numpro = number of boxes being defined (integer)
blow1 = margin added to box to ensure that any desired nodes
lie inside the box, not on the boundary.
(defined in metres) - should be much smaller than grid size.
Lines 2 to (4*numpro)+1:
For each of numpro boxes, there are 4 input lines which describe the
box, including the coordinates of the box starting at the top left
corner and going counter-clockwise, and the material to be assigned
to the box.
Line (i): i_ref_space
i_ref_space - type of definition
0 = absolute X,Z position in metres (joins corners with straight lines)
1 = defined by grid location nx and nz (this will follow any
curvature in grid lines)
Line (ii): xx1 xx2 xx3 xx4
If i_ref_space = 0 then
xx1 = X coordinate of top left corner
xx2 = Z coordinate of top left corner
xx3 = X coordinate of bottom left corner
xx4 = Z coordinate of bottom left corner
If i_ref_space = 1 then
xx1 = nx coordinate of top left corner
xx2 = nz coordinate of top left corner
xx3 = nx coordinate of bottom left corner
xx4 = nz coordinate of bottom left corner
Line (iii): xx5 xx6 xx7 xx8
If i_ref_space = 0 then
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
If i_ref_space = 1 then
xx5 = nx coordinate of bottom right corner
xx6 = nz coordinate of bottom right corner
xx7 = nx coordinate of top right corner
xx8 = nz coordinate of top right corner
Line (iv): ma
ma = number of mechanical property set (as defined in mecprop_i)
to be assigned to this box
Example:
3 10. numpro,blowl (blow for lag. mat. areas)
0
-2.0e+06 2.0e+06 -2.0e+06 -2.0e+06
2.0e+06 -2.0e+06 2.0e+06 2.0e+06
4 crust is b*fspar
0
-1.0e+02 -3.49e+04 -1.0e+02 -5.559e+04
2.999e+05 -5.559e+04 2.699e+05 -3.49e+04
1 mantle
1
-0.5e+00 -27.2e+00 -0.5e+00 -55.2e+00
201.2e+00 -55.2e+00 201.2e+00 -27.2e+00
1 mantle
allocates DEFAULT mechanical sets to eulerian grid 1 'crust or mechanical grid'
L1: number of boxes, tolerance for misfit of lagrangian and eulerian material areas.
==================
File: mecset_L_1_i
Unit: 345
Purpose: Defines the placement of the different materials (which were defined
in mecprop_i) on the Lagrangian grid.
Format/Content:
[NOTE: same format as mecset_E_1_i, except that boxes are defined on the
lagrangian grid.]
Line 1: numpro blowl
numpro = number of boxes being defined (integer)
blow1 = margin added to box to ensure that any desired nodes
lie inside the box, not on the boundary.
(defined in metres) - should be much smaller than grid size.
Lines 2 to (4*numpro)+1:
For each of numpro boxes, there are 4 input lines which describe the
box, including the coordinates of the box starting at the top left
corner and going counter-clockwise, and the material to be assigned
to the box.
Line (i): i_ref_space
i_ref_space - type of definition
0 = absolute X,Z position in metres (joins corners with straight lines)
1 = defined by grid location nx and nz (this will follow any
curvature in grid lines)
Line (ii): xx1 xx2 xx3 xx4
If i_ref_space = 0 then
xx1 = X coordinate of top left corner
xx2 = Z coordinate of top left corner
xx3 = X coordinate of bottom left corner
xx4 = Z coordinate of bottom left corner
If i_ref_space = 1 then
xx1 = nx coordinate of top left corner
xx2 = nz coordinate of top left corner
xx3 = nx coordinate of bottom left corner
xx4 = nz coordinate of bottom left corner
Line (iii): xx5 xx6 xx7 xx8
If i_ref_space = 0 then
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
If i_ref_space = 1 then
xx5 = nx coordinate of bottom right corner
xx6 = nz coordinate of bottom right corner
xx7 = nx coordinate of top right corner
xx8 = nz coordinate of top right corner
Line (iv): ma
ma = number of mechanical property set (as defined in mecprop_i)
to be assigned to this box
Example:
3 10. numpro,blowl (blow for lag. mat. areas)
0
-2.0e+06 2.0e+06 -2.0e+06 -2.0e+06
2.0e+06 -2.0e+06 2.0e+06 2.0e+06
4 crust is b*fspar
0
-1.0e+02 -3.49e+04 -1.0e+02 -5.559e+04
2.999e+05 -5.559e+04 2.699e+05 -3.49e+04
1 mantle
1
-0.5e+00 -27.2e+00 -0.5e+00 -55.2e+00
300.0e+00 -55.2e+00 300.0e+00 -27.2e+00
1 mantle
mechanical set numbers allocated to lgrangian grid. these are regridded
onto eulerian grid 1 when L grid moves.
==================
File: microfem_b_i
Purpose: This is the input file for beam control. It contains information
about beam geometry, strength and loading. (See also Juliet's original
beam notes in Notes 1 book, pages 32-36,134-135)
From the description of microfem input files, by Juliet:
> Beam description - The base of the finite element model is modelled as a
> beam, so that bending (flexure) and isostatic equilibrium (buoyancy -
> balance of forces as though over an enclosed liquid) can be calculated. It
> can either be a single beam or a broken beam. A broken beam will allow a
> discontinuity at the break. The point where the beam is broken is called
> the singularity. In the case of a broken beam, we say there are 2 beams.
> Beam 1 is the LEFT-hand beam (has node numbers > singularity) and is the
> active beam. Beam 2 is the RIGHT-hand beam. Each beam has an end1 at the
> left-hand end, and an end2 at the right-hand end. The variable 'nbreak'
> is the node number of the singularity and is located at end2 of beam2
> (beam nodes are numbered 1 -> # eulerian nodes in horizontal (y)
> direction). The displacement of end2 of beam2 can be slave to that of end1
> of beam1, so any point load must be located >nbreak (usually at nbreak+1).
> Each end of each beam has its own degree-of-freedom (dof) specifications.
> Dof1 is a deflection or force, dof2 is a rotation or moment/torque:
> dof1 = 0 -> imposed deflection
> dof2 = 0 -> imposed rotation (slope)
> dof1 = 1 -> imposed force
> dof2 = 1 -> imposed moment/torque
> A dof=-1 can be specified for dof1 of end2 of beam2, and indicates that
> it is a slave to end1 of beam1 with a ratio rslav.
Format/Content:
Line 1: key
key = title which describes this model run
A70 format
Line 2: nbreak rslav rau1 rau2
Format: i5, 3e9.2
nbreak = node number of end of left-hand beam, where break occurs
rslav = ratio of deflection of slave beam/master beam
rau1 = density of material above beam (crust?)
rau2 = density of material below beam (mantle?)
Line 3: ibbc(1,1) ibbc(1,2) ibbc(1,3) ibbc(1,4)
Format: 4i2
ibbc(1,1) = dof1 specification for end1 of beam 1
ibbc(1,2) = dof2 specification for end1 of beam 1
ibbc(1,3) = dof1 specification for end2 of beam 1
ibbc(1,4) = dof2 specification for end2 of beam 1
dof1 = 0 -> imposed deflection
dof2 = 0 -> imposed rotation (slope)
dof1 = 1 -> imposed force
dof2 = 1 -> imposed moment/torque
See above for more explanation of dof1 and dof2 values.
Line 4: bbc(1,1) bbc(1,2) bbc(1,3) bbc(1,4)
Format: 4e9.2
Values of imposed deflection, rotation, force or moment (per timestep
increments) corresponding to dof's above. The
forces on the beam are cumulative - concentrated loads are specified
below, changing mass distribution and buoyancy are calculated as model
progresses.
Line 5: ibbc(2,1) ibbc(2,2) ibbc(2,3) ibbc(2,4)
Format: 4i2
ibbc(1,1) = dof1 specification for end1 of beam 2
ibbc(1,2) = dof2 specification for end1 of beam 2
ibbc(1,3) = dof1 specification for end2 of beam 2
ibbc(1,4) = dof2 specification for end2 of beam 2
dof1 = 0 -> imposed deflection
dof2 = 0 -> imposed rotation (slope)
dof1 = 1 -> imposed force
dof2 = 1 -> imposed moment/torque
See above for more explanation of dof1 and dof2 values.
Line 6: bbc(2,1) bbc(2,2) bbc(2,3) bbc(2,4)
Format: 4e9.2
Values of imposed deflection, rotation, force or moment (per timestep
increments) corresponding to dof's above.
Line 7: sflex
Format: e9.2
sflex = scaling factor for flexural rigidity (values of flxrig
below will be multiplied by this factor)
Next N lines (Lines 8 to 7+N): flxrig(k),k=1,nodb-1
where N= (number of nodes in beam / 6) rounded up
Format: 6e9.2
This array contains the flexural rigidity of every ELEMENT of the beam,
and will be multiplied by sflex (above).
Next N lines (Lines 8+N to 7+2N): cfb(2*k-1),k=1,nodb
where N = (number of nodes in beam / 6) rounded up
Format: 6e9.2
This array contains the concentrated load (force) at each NODE of the beam.
Next N lines (Lines 8+2N to 7+3N): cfb(2*k),k=1,nodb
where N = (number of nodes in beam / 6) rounded up
Format: 6e9.2
This array contains the concentrated torque at each NODE of the beam.
(not usually .ne.0)
Next line (Line 8+3N): ntfb
Format: i2
ntfb = number of points in the beam time function which follows.
Next lines (Lines 9+3N to end): tb(k),vtfb(k),k=1,ntfb
Format: 2e9.2
beam time function
tb(k) = time in seconds
vtfb(k) = factor to multipy the concentrated load by
To pull the beam down more every
timestep you need to have an increasing time function. The one below
specifies that every timestep the load grows by the value in the
cfb(2k-1) array above.
Example:
m0 b
101 1.00e+00 2.70e+03 3.30e+03
1 1 0 0
0.00E+00 0.00E+00 0.00E+00 0.00E+00
0 0-1 1
0.00E+00 0.00E+00 0.00E+00 0.00E+00
1.00e-01
.10E+24 .10E+24 .10E+24 .10E+24 .10e+24 .10e+24
.10E+24 .10E+24 .10E+24 .10E+24 .10e+24 .10e+24
.10E+24 .10E+24 .10E+24 .10E+24 .10e+24 .10e+24
.10E+24 .10E+24 .10E+24 .10E+24 .10e+24 .10e+24
.10E+24 .10E+24 .10E+24 .10E+24 .10e+24 .10e+24
.10E+24 .10E+24 .10E+24 .10E+24 .10e+24 .10e+24
.10E+24 .10E+24 .10E+24 .10E+24 .10e+24 .10e+24
.10E+24 .10E+24 .10E+24 .10E+24 .10e+24 .10e+24
.10E+24 .10E+24 .10E+24 .10E+24 .10e+24 .10e+24
.10E+24 .10E+24 .10E+24 .10E+24 .10e+24 .10e+24
.10E+24 .10E+24 .10E+24 .10E+24 .10e+24 .10e+24
.10E+24 .10E+24 .10E+24 .10E+24 .10e+24 .10e+24
.10E+24 .10E+24 .10E+24 .10E+24 .10e+24 .10e+24
.10E+24 .10E+24 .10E+24 .10E+24 .10e+24 .10e+24
.10E+24 .10E+24 .10E+24 .10E+24 .10e+24 .10e+24
.10E+24 .10E+24 .10E+24 .10E+24 .10e+24 .10e+24
.10E+24 .10E+24 .10E+24 .10E+24 .10e+24 .10e+24
.10E+24 .10E+24 .10E+24 .10E+24 .10e+24 .10e+24
.10E+24 .10E+24 .10E+24 .10E+24 .10e+24 .10e+24
.10E+24 .10E+24 .10E+24 .10E+24 .10e+24 .10e+24
.10E+24 .10E+24 .10E+24 .10E+24 .10e+24 .10e+24
.10E+24 .10E+24 .10E+24 .10E+24 .10e+24 .10e+24
.10E+24 .10E+24 .10E+24 .10E+24 .10e+24 .10e+24
.10E+24 .10E+24 .10E+24 .10E+24 .10e+24 .10e+24
.10E+24 .10E+24 .10E+24 .10E+24 .10e+24 .10e+24
.10E+24 .10E+24 .10E+24 .10E+24 .10e+24 .10e+24
.10E+24 .10E+24 .10E+24 .10E+24 .10e+24 .10e+24
.10E+24 .10E+24 .10E+24 .10E+24 .10e+24 .10e+24
.10E+24 .10E+24 .10E+24 .10E+24 .10e+24 .10e+24
.10E+24 .10E+24 .10E+24 .10E+24 .10e+24 .10e+24
.10E+24 .10E+24 .10E+24 .10E+24 .10e+24 .10e+24
.10E+24 .10E+24 .10E+24 .10E+24 .10e+24 .10e+24
.10E+24 .10E+24 .10E+24 .10E+24 .10e+24 .10e+24
.10E+24 .10E+24
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00-0.00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00 .00E+00 .00E+00 .00e+00
.00E+00 .00E+00 .00E+00
2
.00E+00 1.00E+00
1.00E+38 1.00E+00
-2.50E+04 2.00E+01
==================
File: microfem_c_i
Unit: 1
Purpose: This is the overall control file for the model run. It determines
how many timesteps will run, the length of each timestep, and which
timesteps will output information for post-processing.
Format/Content:
Line 1: key
key = title which describes this model run
A70 format
Line 2: nstep dtnot grav
Format: free format
nstep = number of timesteps
dt = used to be where length of each timestep specified, not used now
dt is now calculated from vdriv0 and upstep, which are input below.
grav = force due to gravity (absolute value)
Line 3: junk1
Format: (e11.4)
not used
Line 4: itmax er0 ermult
Format: (i5,2e9.2)
not used here, the values in NL_i override these
Line 5: nso ktops
Format: free format
nso = number of timesteps at which to save outputs
ktops = timestep interval at which to save surface topography (we never
plot this so it could probably be increased) - surface topography
is saved in microfem_topo0_s and microfem_topo1_s
NOTE: it's not exactly clear, but I think the code will only write
out the topography when the current timestep is one of nstepo(k)
[see next input line] AND the current timestep is a multiple of
ktops. -BL
Line 6: nstepo(k),k=1,nso
Format: free format
nstepo(k) = Timestep numbers at which to save outputs
Line 7: xsingl xsingr vconv
Format: free format
Juliet says: "These have got something to do with the velocity
singularity (where the velocity along the base goes to 0) moving.
I can't remember the details."
There are two alternative (i.e. mutually exclusive) methods to
specify relationships among convergence velocity (vdrive), length
of timestep (dt) and convergence displacement per timestep
(upstep).
METHOD 1 has fixed upstep (from line 8) but allows velocity to
vary according to a time function.
METHOD 2 has fixed timestep (dtfixed) (see below) but allows
convergence per timestep (upstep) to vary according to a time
function.
All versions of code up to and including TMO4.3 (circa Jan 13
2000) allow only Method 1 and only have corresponding input.
Version TMO4.4 and up (versions created after May 11, 2000) allow
both options.
Line 8: vdriv0 upstep
Format: free format
vdriv0 = basic convergence velocity (m/sec)
upstep = convergence displacement (m)
NOTE: These are parameters for Method 1 (see notes above).
These lines are still required even if Method 2 will be used.
There are two alternative (i.e. mutually exclusive) methods to
specify relationships among convergence velocity (vdrive), length
of timestep (dt) and convergence displacement per timestep
(upstep).
METHOD 1 has fixed upstep (from line 8) but allows velocity to
vary according to a time function.
METHOD 2 has fixed timestep (dtfixed) (see below) but allows
convergence per timestep (upstep) to vary according to a time
function.
All versions of code up to and including TMO4.3 (circa Jan 13
2000) allow only Method 1 and only have corresponding input.
Version TMO4.4 and up (versions created after May 11, 2000) allow
both options.
In Method 1 dt (length of timestep in seconds) is calculated from
dt = upstep/vdrive
where vdrive = vdrive0 X vatime
(vatime specified in Lines 10+)
This system of specifying convergence and dt allows the dt to
change during the model run, so that the convergence velocity
can vary. The amount of convergence (upstep) per timestep
remains constant. This must be used with care because dt can
become very large, or small, depending on vatime and can
strongly modify the accuracy of thermal calculations and erosion
if dt is large.
[Note that true quiescent periods (vatime=0) are also possible
using Method 1. In this case the mechanical f.e. calculations
are made.
notecto = 0 no quiescent period
notecto = 1 quiescent period
see below for corresponding input in file notecto_i]
Line 9: nvetim
Format: free format
nvetim = number of points in velocity time function
Next nvetim lines (Lines 10 to 9+nvetim): atime(k),vatime(k), k=1,nvetim
Format: free format
These points specify a time function controlling the convergence velocity.
The interpolated vatime is multiplied by vdriv0 to calculate the velocity
at a given time. If vatime=0 then no convergence is desired (quiescent
period). In this case, dt = dtquiet (specified below).
Next line (Line 10+nvetim): dtquiet
Format: free format
dtquiet = timestep length in quiescent period (see explanation just above)
The next lines define the parameters for Method 2. Inputs for
Method 2 are not required if Method 1 is being used. Note,
however, there must be an implied 0 (empty line) at the end of the
input file to ensure that nvetimn is read as 0.
Next line: nvetimn
nvetimn = 0 means use Method 1 as specified above
nvetimn > 0 means we are using Method 2 and nvetimn is the
number of points in the convergence time function.
Next nvetimn lines: atimen(k), vatimen(k), k=1,nvetimn
Format: free format
These points specify a time function of the convergence per
timestep. Note that the value of convergence is given,
not a scaling variable.
Next line: dtfixed (only read if nvetimn is not 0)
Format: free format
dtfixed = timestep length for Method 2 calculation
In Method 2 the convergence velocity is calculated from
vatimen/dtfixed. It is therefore possible to vary the
convergence velocity by choosing vatimen(k).
For example, a near quiescent state can be achieved by using
vatimen(k) less than 0.000001 x upstep0 or max(vatimen(k)). In
this case the amount of convergence per timestep and the
convergence velocity can be arbitrarily small to achieve
quiescence. Note that the mechanical (fe) calculation is still
made. *** Convergence criterion will need updating. ****
========================
File: microfem_compile_i
Unit: 710
Purpose: Controls the scripting behaviour of the program ("scripting" means
recording the user's choices in the interactive part of the
program).
Format/Content:
Line 1: i_record_script
Format: free-format integer
i_record_script = 0 means program records interactive input in the
file microfem_compile_s
= 1 means the program will read interactive input
from the file microfem_compile_s
==================
File: microfem_e_i
Purpose: This is the input file for erosion and sedimentation control.
Notes: Juliet says: "if ierost=0, then the amount of erosion
of the surface is calculated using rate, ertime and erspac. This is
somewhat complicated and confusing, and Chris has a summary on the
radiator to the left of adder's console. The summary of it all is:
To INCREASE erosion -> INCREASE erspac
(or) -> DECREASE ertime"
[end of quote from Juliet]
Bonny's addendum: Chris's summary is no longer on the radiator to
the left of Adder's console (office got painted).
Format/Content:
Line 1: key
key = title which describes this model run
A70 format
Line 2: ierost
Format: i2
ierost = 0 -> no erosion of top surface
ierost = 1 -> erosion of top surface proportional to topography
ierost = 2 -> erosion proportional to uplift velocity
ierost = 3 -> erosion proportional to slope
IF ierost = 1 or 2 then the next lines are as follows...
Line 3: rate
Format: e14.7
rate = erosion rate in 1/time
Juliet says: due=dt*h*erspac(now)/(rate*ertime(now))
Line 4: nertim
Format: i2
nertim = number of points in erosion time function
Next nertim lines (Lines 5 to 4+nertim): etime(k),ertime(k) k=1,nertim
Format: 2e14.7
This is the erosion time function, which controls how erosion varies
with time.
etime(k) = time value (units?)
ertime(k) = ??
Next line (Line 5+nertim): nerspa
Format: (i2)
nerspa = number of point in erosion space function
Next nerspa lines: espac(k),erspac(k) k=1,nerspa
Format: (2e14.7)
This is the erosion space function, which controls how erosion varies
across the surface of the model. It should cover the length of the
Eulerian grid.
espac(k) = time value (units?)
erspac(k) = node number (?)
END of section to include if ierost = 1 or 2
IF ierost = 3 then the next lines are as follows...
Line 3: rate
Format: e14.7
Line 4: nertim
Format: i2
Next nertim lines: etime(k), ertime(k)
Format: 2e14.7
Next line: er3_href, er3_alphaplus, er3_alphaminus, er3_stablax
Format: free format
Next line: ner3_x
Format: i2
Next ner3_x lines: er3_x(k), er3_fx(k)
Format: 2e14.7
END of section to include if ierost = 3
Next line: ierosb
Format: (i2)
ierosb = 0 -> no flux of material leaving base
ierosb = 1 -> having material leaving through base
IF iersob = 1 then include the next line:
Next line: nbl nbr
Format: (2i5)
nbl = node number of left-hand side of basal flap
nbr = node number of right-hand side of basal flap
The flap is a region in the base where material can
exit. The velocities in this area have a vertical as
well as a horizontal component, specified in the input
file microfem_g0_i. The endpoints of the flap should have
a z- component = 0, while all the points in between
should have a non-zero z-component.
Example:
m0 e
1 ierost=1 (eros prop to topo) =2 (eros prop to uplift vel)
4.7304000e+11 =rate model dt = 15000*3.1536e7 (attempt at total erosion)
6
0.0000000e+00 1.6666667e+02 etime, ertime
5.5203700e+14 1.6666667e+02 end of model ts 1167
1.5750000e+15 1.6666667e+02 116.7 gives erosion time constant=1.75my
1.7305500e+15 1.6666667e+02
1.7305580e+15 1.6666667e+02
6.0000000e+25 1.6666667e+02
0
WARNING: ierost=1 - erosion=(dt/(rate*ertime))*topo
Make sure that (dt/rate*ertime).le.1
dt may be large during a quiescent period, and ertime must
be adjusted!!
ierost=2 - do not use when model has quiescent period
======================
File: microfem_exhum_i
Purpose: Defines parameters for exhumation map output.
Format/Content:
Line 1: nsavex1 nsavex2 nlefex nrigex timeex
Format: free format
nsavex1 = interval at which to save exhumation maps (in timesteps)
if nsavex1 .le. 0, then exhumation maps will be saved at
the same output timesteps specified in microfem_c_i
nsavex2 = not used
nlefex = element number of left-hand boundary of area to look for
particles (in first layer)
nrigex = element number of right-hand boundary of area to look
for particles (in first layer)
nlefex2 = element number of left-hand boundary of area to look for
particles (in second layer)
nrigex2 = element number of right-hand boundary of area to look
for particles (in second layer)
timeex = time (in seconds) of start of exhumation phase
(maps saved after this time)
Line 2: temp1 temp2 temp3
Format: free format
temp1 = first temperature for cooling history
temp2 = second temperature for cooling history
temp3 = third temperature for cooling history
(The time and depth at which the particles cool past these three
temperatures are saved in the output file microfem_exhum_s.)
Example:
0 0 26 126 0.000000e+00 time (sec) of start of exhumation phase (maps saved)
nsavex1,nsavex2(not used),nlefex,nrigex,timeex
nsavex1 .le. 0 means save exhumation maps at timesteps specified in microfem_c_i
===================
File: microfem_g0_i
Purpose: For the thermo-mechanical code, this file specifies the velocity
boundary conditions for the Eulerian grid. In the original microfem
code, this file also defined the Eulerian grid geometry and the
material properties. The version of this input file for the
thermo-mechanical code still contains lines which define the Eulerian
grid, but they are not used. The Eulerian grid is defined in the
save_grid file instead.
Format/Content:
Line 1: key
Format: A70
key = title - not used by the program, so can be what you want
Line 2: nperel nye nze
Format: (3i4)
nperel = number of nodes per element (restricted to 4)
nye = number of ELEMENTS in horizontal direction
nze = number of ELEMENTS in vertical direction
nye and nze are temporary input names. Throughout the program, the number
of NODES in the horizontal and vertical direction for the Eulerian
grid are commonly stored in igy0 and igz0.
Line 3: scaly scalz
Format: (2e11.4)
scaling factors for the next 4 sets of input (nodal coordinates)
Converts from km to m
Line 4: kcode
Format: a1
kcode = T -> the next few lines specify the location of nodes along
top boundary of Eulerian grid
= B -> the next few lines specify the location of nodes along
bottom boundary of Eulerian grid
= L -> the next few lines specify the location of nodes along
left boundary of Eulerian grid
= R -> the next few lines specify the location of nodes along
right boundary of Eulerian grid
= E -> end of Eulerian grid specification
Next lines: n1 b1 c1 ien1
Format: (i4,2e11.4,i2)
At least two lines of this format are required.
n1 = position of node along the specified boundary (e.g. 1 for first node,
2 for second node, etc.)
b1 = y coordinates of the node (in user units - will be scaled by scaly)
c1 = z coordinates of the node (in user units - will be scaled by scalz)
ien1 = 0 -> this is the end of this block of specifications
= 1 -> more lines of node specifications to follow
Next lines: Repeat lines 4 and following for bottom, left and right
boundaries of grid specification, ending with a line which
contains "E", to end the grid specification.
[The next section specifies the Eulerian grid boundary conditions
using the same kind of format that was used to specify the Eulerian grid
(above).]
Line: kcode
Format: a1
kcode = T -> the next few lines specify the location of nodes for
boundary conditions along the top of the Eulerian grid
= B -> the next few lines specify the location of nodes for
boundary conditions along the bottom of the Eulerian grid
= L -> the next few lines specify the location of nodes for
boundary conditions along the left side of the Eulerian grid
= R -> the next few lines specify the location of nodes for
boundary conditions along the right side of the Eulerian grid
= E -> end of boundary condition specification
Next lines: n1 b1 c1 ien1 iby ibz ien1
Format: (i4,2e11.4,i2)
At least two lines of this format are required.
n1 = position of node along the specified boundary (e.g. 1 for first node,
2 for second node, etc.)
b1 = horizontal velocity in m/ts
c1 = vertical velocity in m/ts
iby = degrees of freedom in the y direction (0=constrained, 1=free)
ibz = degrees of freedom in the z direction (0=constrained, 1=free)
ien1 = 0 -> this is the end of this block of specifications
= 1 -> more lines of node specifications to follow
Next lines: Repeat lines 4 and following for bottom, left and right
boundary conditions, ending with a line which
contains "E" to end.
Next line: iflap
Format: (i2)
iflap - indicates how to treat an opening in the base where mass can exit.
= 0 means there is no mass leaving the base of the grid
= 1 means there is mass leaving the base, and bottom erosion needs
to be specified in microfem_e_i.
Next line: ntflap
Format: i2
ntflap - number of points in flap time function (this value is only
read if nflap .gt. 0)
Next ntflap lines: tflap(k), vtflap(k) k=1,ntflap
Format: (2e9.2)
[These values specify the flap time function and are only read if
iflap.gt.0]
tflap(k) - time (in seconds?)
vtflap(k) - angle of rotation of the basal flap (degrees)? at time tflap(k)
NOTE (from Juliet): If vtflap(k)=0 for all k, the flap stays open and mass
exits governed by the velocities specified in the Eulerian bottom velocity
boundary conditions above. If vtflap varies, then the flap rotates
open/closed, with the amount of rotation specified by vtflap (in degrees?).
More or less mass exits depending on whether the flap is more or less open.
I can't remember exactly how this all works. In some cases the velocity at
the left end of the flap is kept tangent with the slope of the base, and as
the base is pulled down by concentrated loading or increased mass, more
material exits.
Next line: mint
Format: (i5)
mint - number of interfaces to regrid between (usually 2 - top and bottom)
Deflections are calculated at these layers and then interpolated between.
Next mint lines: nnodl0(k) k=1,mint
Format: (10i5)
nnodl0(k) = location of interfaces (row number)
Example:
g
4 200 27
1.0000e+03 1.0000e+03
T
1 0.0000e+00 0.0000e+00 1 |
201 6.0000e+02 0.0000e+00 0 |
B |
1 0.0000e+00-5.5588e+01 1 | NOT USED, Eulerian grid specified
201 6.0000e+02-5.5588e+01 0 | externally from save_grid
L |
1 0.0000e+00 0.0000e+00 1 |
28 0.0000e+00-5.5588e+01 0 |
R |
1 6.0000e+02 0.0000e+00 1 |
28 6.0000e+02-5.5588e+01 0 |
E
B
1-0.0000e+02 0.0000e+00 0 0 1 Basal boundary conditions
101-0.0000e+02 0.0000e+00 0 0 1 displacements in metres/timestep
102-1.5000e+02 0.0000e+00 0 0 1 These ARE used!
201-1.5000e+02 0.0000e+00 0 0 0
L
1-0.0000e+02 0.0000e+00 0 0 1
28-0.0000e+02 0.0000e+00 0 0 0
R
1-1.5000e+02 0.0000e+00 0 0 1
28-1.5000e+02 0.0000e+00 0 0 0
T
1 0.0000e+00 0.0000e+00 1 1 1
201 0.0000e+00 0.0000e+00 1 1 0
E
1
2
.00E+00 .00E+00
.10E+38 .00E+00
2
1 28 0 0 0 0 0 0 0 0
===================
File: microfem_g1_i
Purpose: Defines the initial positions of particles to be traced during the
model run. [NOTE: In the original microfem code, this was the input
file for the Lagrangian grid specifications - grid geometry, particle
information, material properties. For the thermo-mechanical code,
the grid geometry is now specified in save_grid, and the
material properties are specified in mecprop_i. But the particle
positions are still specified in microfem_g1_i.]
Format/Content:
Line 1: blow shift_left_lgrid
Format: (2e9.2)
blow - a height to increase element thickness by (usually some fraction
of element thickness). Sometimes lagrangian nodes will fall outside
Eulerian grid. Normally this means the particles have left the
system, and they're flagged as removed particles (see subroutine
track). Sometimes, however, they are just barely outside the
Eulerian grid, and for numerical (not physical) reasons. 'Blow' lets
you 'blow-up' the top row of Eulerian elements to catch these
particles.
shift_left_lgrid - a distance (in metres) to shift the Lagrangian
grid to the left. (positive value = shift to the left;
negative value = shift to the right). This shift is only done
right after the Lagrangian grid is created. (This gives a
way to extend the Lagrangian grid to the left.)
NOTE: you should make sure that shift_left_lgrid is
equivalent to an integral number of grid cells so that the
Eulerian and Lagrangian grids will match up correctly in
the initial stages of the model.
[Following 2 sets of input lines refer to particles. Particles allow you
to track how rocks which started at a certain location in the grid
would have been displaced and deformed during the model run. You
insert particles by specifying their column and row location in the
Lagrangian grid.]
Line 2: ng2 kt2 dens_ptt
Format: 2i5, e12.4
ng2 - number of particles to track
kt2 - ts interval for saving particle output (is this used?? BL)
dens_ptt - density used to convert depth of burial to lithostatic pressure
Next ng2 lines: icoln,ilin k=1,ng2
Format: (2i5)
icoln - column location of particle
ilin - row location of particle
Next line: junk1 junk2
Format: (2i5)
These values aren't used for anything - don't know why they're still
in the input file.
Next line: ieula1 ieula2 delada xjunk1 xjunk2
Format: (2i3,3e11.4)
ieula1 = ? (not used anywhere in the code)
ieula2 = ? (not used anywhere in the code)
delada = ? (not used anywhere in the code)
xjunk1 = ? (not used anywhere in the code)
xjunk2 = ? (not used anywhere in the code)
Example:
4.00e+02 1.50d+06
30 1 .2700E+04
75 15
75 21
75 27
85 15
85 21
85 27
95 15
95 21
95 27
105 15
105 21
105 27
115 15
115 21
115 27
125 15
125 21
125 27
145 15
145 21
145 27
165 15
165 21
165 27
185 15
185 21
185 27
205 15
205 21
205 27
0 1000
0999 1.0000e+10 1.0000e+10 1.0000e+10
=======================
File: microfem_master_i
Purpose: Determines whether the thermo-mechanical code runs a model or post-
processings an existing models run.
Format/Content:
Line 1: mexec
Format: free format
mexec = 0 -> run
= 1 -> postprocess (plot)
=========================
File: microfem_restartc_i
Purpose: This file controls whether a job is a restart of a previous version
and also specifies how often restart information should be saved etc.
A "restart" of a run will read information from a file called
microfem_restart_s to set its initial internal state.
Format/Content:
Line 1: irestart krestart
Format: (i2 i5)
irestart = 0 -> this run is not a restart (i.e. start job at ts 0)
= 1 -> restart previous job; read information from
microfem_restart_s
krestart - timestep interval for saving restart information
Example:
0 200
=======================
File: microfem_tuning_i
Purpose: This file contains various parameters to tune the model run, dealing
mainly with numerical techniques.
Format/Content:
Line 1: ravog tzer
Format: (2e11.4)
ravog - universal gas constant
tzer - conversion factor: difference between degrees Celsius and
degrees Kelvin
Line 2: iblay imethr icomps izero itherm
Format: (6i2)
iblay - do you want to do horizontal interpolation of grid everywhere or
skip the boundary (basal) layer? (usually do it everywhere)
iblay = 1 forces a constant boundary layer thickness, so keep
it set at 0.
imethr - choice of method or regridding (eg. interpolate only in z, in x,
or in both?)
icomps - ?
izero - ?
itherm - ?
Line 3: icbl icbr ictl ictr
Format: (4i2)
These allow you to 'pin' the corner nodes if you want.
Line 4: irknl
Format: (i2)
irknl = 1 -> use a Runge_Kutta scheme to find velocities.
= 0 -> no
Line 5: iead,icoupling,ieadl
Format: (3i3)
iead - ?
icoupling = 0 -> no coupling between the Lagrangian and Eulerian
mechanical properties
= 1 -> Lagrangian mechanical properties are regridded onto
the Eulerian grid (Note that thermal properties are
regridded automatically)
ieadl - ?
Line 6: ifpa
Format: (i2)
ifpa = 0 -> reuse result_track
= 1 -> make new result_track
> The first time you run the program with given grid inputs, the Lagranian
> and Eulerian grids must be matched. This takes time. When they are
> matched, a file 'result_track' is created which gives the matching.
> Every time you run the program with the same initial grids, you can tell
> the program not to bother matching (ifpa=0) and just read result_track.
> When you change the grids, you need to remember to set ifpa=1 the first
> time. On the SP2, we leave ifpa=1 because the relative time spent in
> this part of the code is small, and it's just one more thing to try and
> remember.
Line 7: viscap
Format: (e9.2)
viscap - limit to maximum viscosity allowed
Example:
8.3144e+00 2.7300e+02 ravog,tzer
0 2 0 0 2 iblay,imethr,icomps,izero,itherm
1 1 0 0 icbl,icbr,ictl,ictr
0 irknl
0 1 0 iead,icoupling,ieadl
1 ifpa
5.00e+17 viscap
icoupling=0 lag mec props not regridded onto eulerian grid
(no coupling between lagrangian and eulerian mec props)
icoupling=1 lag mec props regridded onto eulerian grid
(thermal properties are regridded automatically, so no flag)
====================
File: new_rheologies_i (unit 293)
Purpose:
Format/Content:
Line 1: ddepflag
ddepflag = 1 means read new rheology parameters
Line 2: n_newrheos
Next n_newrheos lines:
ddepparams(i,j), j=1,10
Example:
0 ! 0 if non active 1 if active if =1 DUCTILE parameters in material cards are overwritten BUT plastic parameters stay the same !
1 ! = number of 'colors' (mechanical sets ) then for each of these enter 10 params
2 6d4 1d22 1d20 1d1 1d20 0 0 0 0 !
! 10 parameters are needed to define the rheology of the new material.
! param1 = type of depth-dependence, 1=exponential; 2=linear
! for param1 = 1:
! if z <= param2, then visc = param3
! if z > param2, then visc = param3 X exp (-(Z-param2)/param4)
! if visc <= param5, then visc = param5
! for param1 = 2:
! if z <= param2, then visc = param3
! if z > param2, then visc = param3 + (param4-param3) x (z-param2)/param5
! if visc <= param6, then visc = param6
====================
File: rewind_beamn_i
Purpose: Not sure what this file is for.
Format/Content:
Line 1: ianaly
Format: free format
ianaly = 0 -> ?
ianaly = 1 -> ?? If you look at the code in Xmicrolib.f (subroutine
beam_anal, you'll see
if(ianaly.eq.1)then
write(*,*)'sinusoid test '
... plus more stuff
ianaly = 2 -> ?? Code in Xmicrolib.f (subroutine beam_anal) says:
if(ianaly.eq.2)then
write(*,*)' CONTAMINATION test '
... plus more stuff
=================
File: thermprop_i
Purpose: Defines thermal properties for a set of materials.
Format/Content:
Line 1: mset
mset = number of materials being defined
Lines 2 to (mset*2)+1:
For each material, we have 2 lines in the input file:
thermo(1) thermo(2) thermo(3) thermo(4) thermo(5) thermo(6)
thermo(7) thermo(8) thermo(9) thermo(10) thermo(11) thermo(12)
where
thermo(1) = mantle flux in watts
thermo(2)* = per volume radioactive production in watts/vol
thermo(3)* = depth of radioactive layer (metres)
thermo(4) = thermal conductivity
thermo(5) = density (but only used to calculate diffusivity, not
used in mechanical code)
thermo(6) = Cp (specific heat)
thermo(7) = radioactive heat production per unit mass
(not fully implemented, so set to zero)
thermo(8) = not used
thermo(9) = not used
thermo(10) = not used
thermo(11) = not used
thermo(12) = not used
NOTE:
The pre-processor does NOT use any of these properties. Also, any values
used in the pre-processor are not carried over to the rest of the model run.
* These parameters are not used anymore; they were used to define the
radioactive heat-producing layer, but now replaced by the boxes defined in
therset_L_1_i and therset_L_1_radio_i.
Thermo(1,1) specifies the mantle heat flux and MUST be set correctly. This
is not specified by the pre-processing phase. Check compatible with
pre-processor mantle heat flux.
Example:
6
3.00e-02 0.00e-06 2.00e+04 2.25e+00 3.00e+03 7.50e+02
0.00e-06 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
3.00e-02 0.00e-06 2.00e+04 2.25e+00 3.00e+03 7.50e+02
0.00e-06 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
3.00e-02 0.00e-06 2.00e+04 2.25e+00 3.00e+03 7.50e+02
0.00e-06 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
3.00e-02 0.00e-06 2.00e+04 2.25e+00 3.00e+03 7.50e+02
0.00e-06 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
3.00e-02 0.00e-06 2.00e+04 2.25e+00 3.00e+03 7.50e+02
0.00e-06 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
3.00e-02 0.00e-06 2.00e+04 2.25e+00 3.00e+03 7.50e+02
0.00e-06 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00
thermal properties (here 12 but you can modify that if wanted)
mantle flux in watts , per volume radioactive production in watts/vol,
depth of radoactive layer in meters, conductivity(NOT difusivity) ,'density'
(for diffusivity computation), Cp, again parameter 2.
===================
File: therset_E_2_i
Unit: 346
Purpose: Defines the placement of the thermal materials (which were defined
in thermprop_i) on the thermal grid.
Format/Content:
[NOTE: This format is just like the format for mecset_E_1_i]
Line 1: numpro blowl
numpro = number of boxes being defined (integer)
[The boxes may be radiogenic as defined in therset_E_2_radio_i]
blow1 = margin added to box to ensure that any desired nodes
lie inside the box, not on the boundary.
(defined in metres) - should be much smaller than grid size.
Lines 2 to (4*numpro)+1:
For each of numpro boxes, there are 4 input lines which describe the
box, including the coordinates of the box starting at the top left
corner and going counter-clockwise, and the material to be assigned
to the box.
Line (i): i_ref_space
i_ref_space - type of definition
0 = absolute X,Z position in metres (joins corners with straight lines)
1 = defined by grid location nx and nz (this will follow any
curvature in grid lines)
Line (ii): xx1 xx2 xx3 xx4
If i_ref_space = 0 then
xx1 = X coordinate of top left corner
xx2 = Z coordinate of top left corner
xx3 = X coordinate of bottom left corner
xx4 = Z coordinate of bottom left corner
If i_ref_space = 1 then
xx1 = nx coordinate of top left corner
xx2 = nz coordinate of top left corner
xx3 = nx coordinate of bottom left corner
xx4 = nz coordinate of bottom left corner
Line (iii): xx5 xx6 xx7 xx8
If i_ref_space = 0 then
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
If i_ref_space = 1 then
xx5 = nx coordinate of bottom right corner
xx6 = nz coordinate of bottom right corner
xx7 = nx coordinate of top right corner
xx8 = nz coordinate of top right corner
Line (iv): ma
ma = number of thermal property set (as defined in thermprop_i)
to be assigned to this box
NOTE that corresponding radioactive properties are specified in the same
order in therset_E_2_radio_i; i.e. the first box in therset_E_2_i has
the radioactive properties of the first set in therset_E_2_radio_i.
Example:
6 10. numpro,blowl (blow for lag. mat. areas)
0
-2.0e+06 2.0e+06 -2.0e+06 -2.0e+06
2.0e+06 -2.0e+06 2.0e+06 2.0e+06
6
0
0.0e+00 0.0e+00 0.0e+00 -2.0e+04
2.00e+06 -2.0e+04 2.00e+06 0.0e+00
1
0
5.00e+05 -5.00e+04 5.00e+05 -5.10e+04
5.01e+05 -5.10e+04 5.01e+05 -5.00e+04
3
0
2.0e+05 -2.1e+04 2.2e+05 -2.8e+04
3.00e+05 -2.6e+04 3.00e+05 -2.4e+04
2
0
-1.0e+02 -3.5e+04 -1.0e+02 -1.8e+06
6.0e+05 -1.8e+06 6.0e+05 -3.5e+04
4
1
6.25e+01 -2.75e+01 6.25e+01 -4.25e+01
1.015e+02 -4.25e+01 1.015e+02 -2.75e+01
5
where are thermal colors in 'thermal grid'
=========================
File: therset_E_2_radio_i
Unit: 746
Purpose: Defines the radiogenic heat production for the areas defined in
therset_E_2_i.
Format/Content:
Line 1: numpro blowl
numpro = number of boxes being defined (integer)
[This number must be the same as numpro in therset_E_2_i]
blow1 = margin added to box to ensure that any desired nodes
lie inside the box, not on the boundary.
(defined in metres) - should be much smaller than grid size.
Next lines:
For each of numpro boxes, define the radiogenic heat production at
the corners of the box, starting from top left and going counterclockwise,
and define a time function for the heat production.
Line (i): ff1 ff2 ff3 ff4
Format: free format
ff1 = radiogenic heat production at top left corner of box
ff2 = radiogenic heat production at bottom left corner of box
ff3 = radiogenic heat production at bottom right corner of box
ff4 = radiogenic heat production at top right corner of box
[NOTE: Inside the box, the value is given by linear interpolation
of the corner values.]
Line (ii): nth_time
Format: free format
nth_time = number of values in the time function
Next nth_time lines: th_mat_time, th_mat_val
Format: free format
[This is the time function for the heat production in this box]
th_mat_time - time (in seconds)
th_mat_valu - strength multiplier (uniformly scales values of
radiogenic heat production in the box)
NOTE that each set corresponds to the equivalent thermal box in therset_E_2_i;
i.e. first box in therset_E_2_i has the radioactive properties of the first set
in therset_E_2_radio_i.
Example:
6 10. numpro,blowl (blow for lag. mat. areas)
0.0e-06 0.0e+00 0.0e+00 0.0e-06
6
0.0000000e+00 1.0000000e+00
5.5203700e+14 1.0000000e+00 end of model ts 1167
1.5750000e+15 1.0000000e+00
1.7305500e+15 1.0000000e+00
1.7305580e+15 1.0000000e+00
6.0000000e+25 1.0000000e+00
2.0e-06 2.0e-06 2.0e-06 2.0e-06
6
0.0000000e+00 1.0000000e+00
5.5203700e+14 1.0000000e+00 end of model ts 1167
1.5750000e+15 1.0000000e+00
1.7305500e+15 1.0000000e+00
1.7305580e+15 1.0000000e+00
6.0000000e+25 1.0000000e+00
0.0e-06 0.0e-06 0.0e-06 0.0e-06
6
0.0000000e+00 1.0000000e+00
5.5203700e+14 1.0000000e+00 end of model ts 1167
1.5750000e+15 1.0000000e+00
1.7305500e+15 1.0000000e+00
1.7305580e+15 1.0000000e+00
6.0000000e+25 1.0000000e+00
0.0e-06 0.0e+00 0.0e+00 0.0e-06
6
0.0000000e+00 1.0000000e+00
5.5203700e+14 1.0000000e+00 end of model ts 1167
1.5750000e+15 1.0000000e+00
1.7305500e+15 1.0000000e+00
1.7305580e+15 1.0000000e+00
6.0000000e+25 1.0000000e+00
0.0e-06 0.0e+00 0.0e+00 0.0e-06
6
0.0000000e+00 1.0000000e+00
5.5203700e+14 1.0000000e+00 end of model ts 1167
1.5750000e+15 1.0000000e+00
1.7305500e+15 1.0000000e+00
1.7305580e+15 1.0000000e+00
6.0000000e+25 1.0000000e+00
0.0e-06 0.0e-06 0.0e-06 0.0e-06
6
0.0000000e+00 1.0000000e+00
5.5203700e+14 1.0000000e+00 end of model ts 1167
1.5750000e+15 1.0000000e+00
1.7305500e+15 1.0000000e+00
1.7305580e+15 1.0000000e+00
6.0000000e+25 1.0000000e+00
wher are thermal colors in 'tghermal grid'
===================
File: therset_L_1_i
Unit: 347
Purpose: Defines the thermal property boxes on the Lagrangian grid.
NOTE: Lagrangian thermal properties will overwrite Eulerian thermal properties
on the part of the thermal grid where the Lagrangian overlaps the thermal
grid. This means that in the region covered by the Lagrangian grid,
the thermal properties will be advected with the mechanical flow.
Elsewhere on the thermal grid, thermal properties are defined by
the Eulerian thermal sets, and these will not be advected.
Format/Content:
[NOTE: This format is just like the format for therset_E_2_i]
Line 1: numpro blowl
numpro = number of boxes being defined (integer)
[The boxes may be radiogenic as defined in therset_L_1_radio_i]
blow1 = margin added to box to ensure that any desired nodes
lie inside the box, not on the boundary.
(defined in metres) - should be much smaller than grid size.
Lines 2 to (4*numpro)+1:
For each of numpro boxes, there are 4 input lines which describe the
box, including the coordinates of the box starting at the top left
corner and going counter-clockwise, and the material to be assigned
to the box.
Line (i): i_ref_space
i_ref_space - type of definition
0 = absolute X,Z position in metres (joins corners with straight lines)
1 = defined by grid location nx and nz (this will follow any
curvature in grid lines)
Line (ii): xx1 xx2 xx3 xx4
If i_ref_space = 0 then
xx1 = X coordinate of top left corner
xx2 = Z coordinate of top left corner
xx3 = X coordinate of bottom left corner
xx4 = Z coordinate of bottom left corner
If i_ref_space = 1 then
xx1 = nx coordinate of top left corner
xx2 = nz coordinate of top left corner
xx3 = nx coordinate of bottom left corner
xx4 = nz coordinate of bottom left corner
Line (iii): xx5 xx6 xx7 xx8
If i_ref_space = 0 then
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
If i_ref_space = 1 then
xx5 = nx coordinate of bottom right corner
xx6 = nz coordinate of bottom right corner
xx7 = nx coordinate of top right corner
xx8 = nz coordinate of top right corner
Line (iv): ma
ma = number of thermal property set (as defined in thermprop_i)
to be assigned to this box
NOTE that corresponding radioactive properties are specified in the same
order in therset_L_1_radio_i; i.e. the first box in therset_L_1_i has
the radioactive properties of the first set in therset_L_1_radio_i.
Example:
4 10. numpro,blowl (blow for lag. mat. areas)
0
-2.0e+06 2.0e+06 -2.0e+06 -2.0e+06
2.0e+06 -2.0e+06 2.0e+06 2.0e+06
6
0
0.0e+00 0.0e+00 0.0e+00 -2.0e+04
2.00e+06 -2.0e+04 2.00e+06 0.0e+00
1
0
5.00e+05 -5.00e+04 5.00e+05 -5.10e+04
5.01e+05 -5.10e+04 5.01e+05 -5.00e+04
3
0
2.0e+05 -2.1e+04 2.2e+05 -2.8e+04
3.00e+05 -2.6e+04 3.00e+05 -2.4e+04
2
where are thermal 'colors' on lagrangian grid. note these are adevcted.
and regridded to 'thermal grid'
=========================
File: therset_L_1_radio_i
Unit: 747
Purpose: Defines the radiogenic heat production of the boxes which were defined
in therset_L_1_i.
Format/Content:
Line 1: numpro blowl
numpro = number of boxes being defined (integer)
[This number must be the same as numpro in therset_L_1_i]
blow1 = margin added to box to ensure that any desired nodes
lie inside the box, not on the boundary.
(defined in metres) - should be much smaller than grid size.
Next lines:
For each of numpro boxes, define the radiogenic heat production at
the corners of the box, starting from top left and going counterclockwise,
and define a time function for the heat production.
Line (i): ff1 ff2 ff3 ff4
Format: free format
ff1 = radiogenic heat production at top left corner of box
ff2 = radiogenic heat production at bottom left corner of box
ff3 = radiogenic heat production at bottom right corner of box
ff4 = radiogenic heat production at top right corner of box
[NOTE: Inside the box, the value is given by linear interpolation
of the corner values.]
Line (ii): nth_time
Format: free format
nth_time = number of values in the time function
Next nth_time lines: th_mat_time, th_mat_val
Format: free format
[This is the time function for the heat production in this box]
th_mat_time - time (in seconds)
th_mat_valu - strength multiplier (uniformly scales values of
radiogenic heat production in the box)
NOTE that each set corresponds to the equivalent thermal box in therset_L_1_i;
i.e. first box in therset_L_1_i has the radioactive properties of the first set
in therset_L_1_radio_i.
Example:
4 10. numpro,blowl (blow for lag. mat. areas)
0.0e-06 0.0e+00 0.0e+00 0.0e-06
6
0.0000000e+00 1.0000000e+00
5.5203700e+14 1.0000000e+00 end of model ts 1167
1.5750000e+15 1.0000000e+00
1.7305500e+15 1.0000000e+00
1.7305580e+15 1.0000000e+00
6.0000000e+25 1.0000000e+00
2.0e-06 2.0e-06 2.0e-06 2.0e-06
6
0.0000000e+00 1.0000000e+00
5.5203700e+14 1.0000000e+00 end of model ts 1167
1.5750000e+15 1.0000000e+00
1.7305500e+15 1.0000000e+00
1.7305580e+15 1.0000000e+00
6.0000000e+25 1.0000000e+00
0.0e-06 0.0e-06 0.0e-06 0.0e-06
6
0.0000000e+00 1.0000000e+00
5.5203700e+14 1.0000000e+00 end of model ts 1167
1.5750000e+15 1.0000000e+00
1.7305500e+15 1.0000000e+00
1.7305580e+15 1.0000000e+00
6.0000000e+25 1.0000000e+00
0.0e-06 0.0e+00 0.0e+00 0.0e-06
6
0.0000000e+00 1.0000000e+00
5.5203700e+14 1.0000000e+00 end of model ts 1167
1.5750000e+15 1.0000000e+00
1.7305500e+15 1.0000000e+00
1.7305580e+15 1.0000000e+00
6.0000000e+25 1.0000000e+00
where are thermal 'colors' on lagrangian grid. note these are adevcted.
and regridded to 'thermal grid'
===============
File: save_grid
Purpose: This file defines the thermal grid (and Eulerian grid? - ask Philippe)
It's not meant to be created by hand. The save_grid file that we've
been using for most of the thermo-mechanical model runs was created by
the thermo-mechanical code (see directions in another document).
Eventually I will create a stand-alone program to generate the
save_grid file.
Example (first few lines only):
grid is
topology 201 55
total number of nodes 11055
1 .0 .0
2 .0 -2058.81485502604
3 .0 -4117.62971005208
4 .0 -6176.44420317708
5 .0 -8235.25942010417
6 .0 -10294.0739132292
7 .0 -12352.8884063542
... etc...