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...