INPUTS FOR SOPALE1 (Serial OPALE) dmtm version [Deep Model Thermal Mechanical with erosion] (code from September 19, 2002) --------------------------------------------------------------------- File: SOPALE1_i read(8,*)template_i template_i = number of lines to skip used for reading header= model templates listed at the beginning of the file do i=1,template_i read(8,'(a1)')modelkey enddo This part skips over the lines which describe the model templates read(8,*)nmodels nmodels = number of models read(8,*)(modelsn(i),i=1,nmodels) list of model numbers to run For each model, we have the following input: read(8,'(a1)')modelkey modelkey = one-character identifier for the model read(8,*)amodel, benchmark model number for this model, benchmark number read(8,*)wres,rres wres = number of timesteps between writing of 2 restart files rres = determines whether to read restart file rres = 0 - don't read restart file rres = 1 - read restart file read(8,*)nsave1 number of saves for Eulerian grid (Maximum number of saves is 200) read(8,*)(isave_1(i),i=1,nsave1) [only needed if nsave1 > 0 ] timesteps at which to save the Eulerian grid read(8,*)nsave2 number of saves for Lagrangian grid (Maximum number of saves is 200) read(8,*)(isave_2(i),i=1,nsave2) [only needed if nsave2 > 0 ] timesteps at which to save the Lagrangian grid read(8,'(a14)')nameout1 14-character prefix for the name of the output files NOTE: there must be NO MORE and NO LESS than 14 characters, and the last two characters must be "/g". Bonny's note: Philippe was assuming that the prefix would always be './Sopale/rXX/g' where XX is 01, 02, 03, etc. Because of the way Philippe wrote the code, you cannot use a prefix of more than or less than 14 characters. Using blanks before or after a shorter prefix will give problems with embedded blanks in filenames. This means that if you want to keep the input files and the output file in the same directory, you would have to use the prefix '././././././/g'. read(8,*)time_nstep,dt,outunt time_nstep = number of timesteps in this model run dt = timestep length outunt = Fortran unit number for output, e.g. 60 read(8,*)bigu,bigt bigu = ~1d7 multipliers for penalizing boundary velocity bigt = ~1d7 multiplier for penalizing the temperature boundary conditions read(8,*)length,height,nx1,ny1,betaxg,betayg,length2,height2,nx2,ny2 length = length of Eulerian rectangle height = height of Eulerian rectangle nx1 = number of rows of Eulerian grid ny1 = number of columns of Eulerian grid betaxg = 1 for regular grid (heterogeneous grids are other values) betayg = 1 for regular grid (heterogeneous grids are other values) length2 = length of Lagrangian rectangle height2 = height of Lagrangian rectangle nx2 = number of rows of Lagrangian grid ny2 = number of columns of Lagrangian grid NOTE: surface is at "height", not 0 nodes are numbered starting at top surface x=0 is at left edge read(8,*)npert npert = number of interfaces to perturb for instability analysis = 0 for sandbox analysis read(8,*)(ipert(i),i=1,npert) [only required if npert > 0] read(8,*)(lpert(i),i=1,npert) [only required if npert > 0] read(8,*)(lpert0(i),i=1,npert) [only required if npert > 0] read(8,*)(apert(i),i=1,npert) [only required if npert > 0] read(8,*)nfpert nfpert = perturbation of forces = 0 for our models read(8,*)afpert,flpert0,flpert,fipert,xyfpert [only required if nfpert > 0] read(8,*)iformul read(8,*)(hists(i),i=1,11) read(8,*)ang_gravity,gravity,pcohes,visminmax(1),visminmax(2),psmooth, nelvismin1, nelvismin2, nelvismin3, visminscale ang_gravity = angle of gravity in degrees The convention is that ang_gravity is positive when measured anti-clockwise from the positive X direction. Gravity directed in positive X direction has ang_gravity=0. gravity = value of gravity assumed to act toward decreasing Z pcohes = lower limit for pressure (i.e. any pressure below pcohes will be set to pcohes) visminmax(1) = minimum viscosity limit (any viscosities below this are set to visminmax(1)) visminmax(2) = maximum viscosity limit (any viscosities above this are set to visminmax(2)) psmooth = 1 for pressure smoothing = anything else, no pressure smoothing nelvismin1 = element number for first region of viscosity modification. Elements from 1 to nelvismin1 use a minimum viscosity of visminmax(1). nelvismin2 = element number for second region of viscosity modification. Elements from 1+nelvismin1 to nelvismin2 use a minimum viscosity of visminmax(1)/visminscale. nelvismin3 = element number for third region of viscosity modification. Elements from 1+nelvismin2 to nelvismin3 use a minimum viscosity of visminmax(1)/visminscale**2. visminscale = scaling factor for visminmax(1) [see descriptions of nelvismin1, nelvismin2, nelvismin3]. NOTE 1: if nelvismin1, nelvismin2, nelvismin3, and visminscale are not specified (as will happen with earlier versions of SOPALE1_i), the code will set them to zero and the minimum viscosity will not be used. To apply the minimum viscosity limit [visminmax(1)] as in previous versions of SOPALE, set nelvismin1, nelvismin2, and nelvismin3 greater than the largest element number in your model. NOTE 2: One of the uses of the minimum viscosity limit is to prevent slumping in a model. However, a disadvantage of this method is that the entire model is affected, when the slumping problem may be only in the surface. The model is even more sensitive to the effect of vismin when using fluid pressure (lambda_pore_prs). Introducing nelvismin1, nelvismin2, nelvismin3, and visminscale was an attempt to stabilize the surface with a minimum viscosity while making a gradual transition to unmodified viscosity further down in the model. So far [as of 2002/09/19] we (CB & BL) haven't had much success with this technique. read(8,*)color_number color_number = number of materials whose mechanical properties are to be defined read(8,*)mixhyp,refsrate,redphi,nredphi plas_fluid_method mixhyp = mixing hypothesis for viscous material properties (viscous rheologies) = 1 viscous rheologies act in serial NOT CODED = 2 viscous rheologies act in parallel. Each viscous rheology is assumed to have the same second invariant of the strain rate (they act in parallel). Stresses are determined from their respective effective viscosities X strain rate. Total second invariant of stress is sum of individual stresses. The total viscous stress is compared with the plastic yield stress to determine whether flow is viscous or plastic. = 3 minimum dissipation Viscous rheologies combine in a way that minimizes dissipation NOT CODED Use mixhyp = 2 refsrate = reference strain rate to start viscosity computation for linear creep NOTE: special method for setting mechanical material properties, used primarily to create a bottom boundary layer. This sets the material properties for the mechanical calculations only, it does NOT change the actual color (material) number of the Eulerian elements. Set nredphi > total number of Eulerian elements to suppress this method. redphi = mechanical material number for all elements greater than nredphi nredphi = element number. Eulerian elements greater that this number will use the properties of material redphi for mechanical calculations. Note: elements are numbered by row, so this can be used to set a basal boundary layer. The mechanical material properties overrides those specified through the normal input process (Quick and dirty method for setting all elements to material redphi - use redphi= 0) CB Comments: redphi and nredphi are used to define special conditions for the mechanical Eulerian element materials. If the element number is greater than nredphi, it will be assigned the mechanical material properties of material number redphi for mechanical calculations. Note that the actual color number (material number) of the Eulerian element is not overwritten, so you will not see the boundary layer material redphi in plots of material colors. Use of redphi and nredphi overrides other methods used to assign mechanical properties. [Example of use is to define a basal boundary layer in sandbox type experiments such that the boundary layer is redefined every timestep according to the Eulerian grid and therefore in a manner that overrides material advection and regridding from the Lagrangian grid.] nredphi should be used with care, remember that it assigns one mechanical material property (that of redphi) to the Eulerian elements. The effect of strain dependent material properties should also be considered. To avoid resetting element properties in this way, set nredphi.gt.max number of Eulerian mechanical elements. plas_fluid_method = method for calculating plastic yielding for a frictional plastic incompressible material using a Drucker-Prager yield, for plane-strain and including various options for the effects of pore fluid pressure (see notes below). plas_fluid_method not specified in input or plas_fluid_method = 0 Defaults to old behaviour in which phi values (see below) are interpreted as effective phi's modified by pore fluid pressure. In subroutine viscosities1, yielding is given by splas = epressure * dsin(phi) + coh plas_fluid_method = 1 Assume lambda_pore_prs = fluid pressure/solid pressure (mean stress) splas = epressure*(1-lambda_pore_prs)*dsin(phi) + coh (small error in last term) plas_fluid_method = 2 Assume lambda_pore_prs = fluid pressure/solid pressure (mean stress) splas = epressure*(1-lambda_pore_prs)*dsin(phi) + coh*dcos(phi) (small error corrected) plas_fluid_method = 3 (see warning below) Assume lambda_pore_prs = fluid pressure/overburden weight (in code, overburden weight is called plitho) splas = (epressure - lambda_pore_prs*plitho)*dsin(phi) + coh (small error in last term) plas_fluid_method = 4 (see warning below) Assume lambda_pore_prs = fluid pressure/overburden weight (in code, overburden weight is called plitho) splas = (epressure - lambda_pore_prs*plitho)*dsin(phi) + coh*dcos(phi) (small error corrected) plas_fluid_method = 5 The effect of hydrostatic fluid pressures on epressure is taken into account. Method assumes that the fluid pressure is equal to the hydrostatic fluid pressure everywhere. Because fluid pressure is hydrostatic, lambda_pore_prs defaults to 1. (in code, hydrostatic weight is called phydro) splas = (epressure - lambda_pore_prs*phydro)*dsin(phi) + coh*dcos(phi) When you use plas_fluid_method=5, lambda_pore_prs1 and lambda_pore_prs2 must be set to 1 in SOPALE1_i. NOTE: The hydrostatic pressure (phydro) is applied to ALL materials in the model. User beware! NOTE: Methods 1 and 3 have errors in the cohesion term. These methods are available for comparison with models that used the old behaviour. Methods 2 and 4 have the correct cohesion term. WARNING! As of August 2002, models using plas_fluid_method = 3 and plas_fluid_method = 4 will not converge because epressure can be less than lambda_pore_prs*plitho, resulting in a negative splas. A fix is still under consideration. ** UPDATE ** [BL 2002/09/18] Using lambda_blend1 and lambda_blend2 to gradually introduce the lambda_pore_prs seems to improve the chances of the model getting past the first timestep. Also, values of lambda_pore_prs greater than 0.4 seem to cause problems. For each material (i=1,color_number), we have the following input: read(8,*)density(i) density(i) = density of this material read(8,*)bulkv(i) bulkv(i) = bulk viscosity NOTE: Plasticity yield criterion updated July 23, 2002 to include fluid pore pressure effects explicitly. Plastic yielding uses the Drucker-Prager criterion 1/2 (J ') = p(1-lambda) X sin(phi) + c 2 where J ' is the second invariant of the deviatoric stress 2 p is the pressure in the solid phi is the internal angle of friction c is the cohesion lambda is the pore pressure ratio (p / p ) fluid solid In subroutine viscosities1, this appears as splas = epressure*(1.0-lambda_pore_prs)*dsin(phi) + coh (Old version was splas=epressure*dsin(phi) + coh) phi, lambda_pore_prs, and coh are material properties specified in input for each material (see below, description of next two lines of input). phi, lambda_pore_prs and coh can be specified to vary with the second invariant of the deviatoric strain, here termed "strain". read(8,*)phi1(i),phi2(i),phie1(i),phie2(i), lambda_pore_prs1(i), lambda_pore_prs2(i), lambda_pore_prse1(i), lambda_pore_prse2(i) phi1(i) = first angle of friction in degrees phi2(i) = second angle of friction in degrees phie1(i) = strain level 1 phie2(i) = strain level 2 NOTE: must have phie1(i) < phie2(i) Diagram: ________________ phi2 /| / | / | / | phi1 ______________/ phie2 (strain level 2) | phie1 (strain level1) lambda_pore_prs1(i) = first pore pressure ratio lambda_pore_prs2(i) = second pore pressure ratio lambda_pore_prse1(i) = strain level 1 lambda_pore_prse2(i) = strain level 2 NOTE: must have lambda_pore_prse1(i) < lambda_pore_prse2(i) ________________ lambda_pore_prs2 /| / | / | / | lambda_pore_prs1 __________/ lambda_pore_prse2 | lambda_pore_prse1 Use of phi and lambda_pore_prs Previous use varied phi to take account of lambda_pore_prs implicitly, i.e. by choosing phi_effective such that sin(phi_effective) = (1.0-lambda_pore_prs)*sin(phi) This method can still be used in which case all inputs for labmda_pore_prs should be set to 0 and the input values of phi are interpreted as effective internal angles of friction. Alternatively, the input phi's are true, material properties, e.g. compatible with Byerlee's Law (30-40 degrees) and the labmda_pore_prs values are effective pore pressure ratios associated with each material. Note fluid pressures are not calculated separately/dynamically in this verison. Lambda_pre_prs is specified as a material property and should be chose to reflect the fluid pressure environment envisioned for this material. lambda_pore_prs = 0.0 dry rock environment, no fluid pressure 0.3-0.4 hydrostatic fluid pressure - depends on the density distribution of the model material close to 1.0 (but not equal to 1.0) lithostatic conditions, strongly overpressured conditions Note that the theory used here is the most simple form. It hids the underlying problem of the best measure of the solid pressure in lambda_pore_prs. The solid pressure may not be lithostatic. This is likely the case in frictional-plastic materials (see Petrini and Podladchikov, J.Met.Geol.v.18, p.67-77,2000). NOTES CONCERNING CHOICE OF FRICTIONAL PLASTIC FAILURE CRITERIA IN SOPALE- METHODS 1 AND 2 FOR INCLUDING PORE FLUI PRESSURE EFFECTS USING THE HUBBERT-RUBEY PORE FLUID PRESSURE RATIO Derivation is 2D and then generalised to Plane Strain Frictional Plastic Failure Criteria Including Effects of Pore Fluid Pressure Note in this section the names of variables/parameters are not the same as those used elsewhere in the guide. The approx conversion is phi = f lambda_pore_prs = L coh = C stress = S (various components) 1) Definition of Hubbert-Rubey Pore Fluid Pressure Ratio As Glen points out Hubbert and Rubey have defined L (or lambda) as L1 = Pf / Sn and as L2 = Pf / Sz where Pf is the pore fluid pressure, Sn is the normal stress component acrossa potential failure plane and Sz is the weight of the overburden ( equivalent to the lithostatic weight rho x g x z in a uniform medium). Here we use the second definition in the same manner as Davis et al 1983. They define the vertical normal traction to be solely the effect of the lithostatic overburden and include a hydrostatic component for submarine problems. To avoid confusion it is probably best to refer to Sz as the weight of the overburden. The confusing term is the 'lithostatic pressure'...which is sometimes used to mean the weight of the overburden and sometimes used to mean the solid pressure. The stress convention is that compressive stresses are positive ( this is not the normal convention...but it is also used by Davis et al and Jaeger and Cook.) 2) Definition of Effective Stress Se(1) = S(1) - Pf, Se(3) = S(3) - Pf , where Se(1), S(1) are the maximum effective principal stress and the maximum principal stress and the Se(3) and S(3) are the corresponding minimum principal stresses. 3) The Coulomb Criterion (Jaeger and Cook, page 95) Shear failure in rock occurs when Mod( T ) = C + S x Tan (f) (1 where S and T are the normal and shear stresses across the failure plane, C is a constant, the cohesion and f is the internal angle of friction ( normally called phi). This can also be expressed in terms of the Mohr circle representation. 1/2 x (S(1)-S(3)) = 1/2 x(S(1)+S(3))x Sin(f) + C x Cos(f) (2 4) Effect of Fluid Pressure on the Coulomb Failure Criterion The Mohr envelope for the Coulomb criterion is modified by pore fluid pressure to Mod (T) = (S - Pf) x Tan (f) + C or Mod (T) = Se x Tan (f) + C where Se is the effective normal stress on the failure plane. Similarly the Mohr circle representation can be written in terms of effective stresses 1/2 x( Se(1)-Se(3)) = 1/2 x (Se(1)+Se(3))x Sin(f) + C x Cos(f) (3 5) Mean Solid Pressure the mean pressure in the solid is P = 1/2 x (S(1)+S(3)) Therefore P = 1/2 x (Se(1)+Se(3)) + Pf (4 Substituting for Se(1)+Se(3) from 4 into 3 1/2 x(Se(1)-Se(3)) = (P - Pf) x Sin(f) + C x Cos(f) (5 Because Se(1)-Se(3) = S(1)-S(3) 5 can be written 1/2 x(S(1)-S(3)) = (P - L x Sz) x Sin(f) + C x Cos(f) This is the two dimensional form of the Coulomb failure criterion taking pore fluid pressure into account through the Hubbert-Rubey pore fluid pressure ratio. 6) Failure Criterion in Plane Strain in Terms of the Stress Invariants If it is assumed that in plane strain that the intermediate principal stress acts in the direction normal to the plane containing the other principal stresses and has a value approximately equal to the solid pressure then both the second invariant of the stress and the second invariant of the deviatoric stress both reduce to Sqrt (I2) = Sqrt (J2D) = 1/2 x (S(1)-S(3)) where Sqrt is the square root. ASIDE This is because for S(2) = P it follows that S(1)-S(3) = 2 x(S(2)-S(3)) = 2 x (S(1)-S(2)) and therefore Sqrt (I2)= 1/6 x [(S(1)-S(3))^2 + (S(3)-S(1))^2 + (S(1)-S(2))^2] becomes = 1/6x[ (1/4 + 1 + 1/4)(S(1)-S(3))] = 1/2 x(S(1)-S(3))....as above Note this is only correct if the intermediate principal stress is out of the plane and equal to the pressure. It neeeds to be shown that this is also true if yielding is written in terms of the second invariant of the deviatoric stress. 7) 'Mohr-Coulomb' Failure Law Montesi and Zuber give the M-C failure law in terms of the invariants as I2 = I1 x Sin(f) + C x Cos(f) their eqn 18 This is close to the form given above but differs in sqrt(I2) vs I2 and I1 vs P. The eqn is not dimensionally correct with I2 .They may choose definitions of I2 as sqrt(I2) and I1 = I1/3 = P END ASIDE Under these circumstances the failure criterion can be written as Sqrt (J2) = Sqrt (J2D) = (P - L x Sz) x Sin(f) + C x Cos(f) where P = 1/3 (I1) = 1/3 x ((S(1) + S(2) + S(3)) The above is the failure criterion when Plas_Fail_ Method = 2 and appears in subroutine viscosities as splas = (epressure - lambda_pore_prs x plitho) x dsin(phi) + coh Note the small error that coh should be coh x dcos(phi) END NOTES read(8,*)coh1(i),coh2(i),cohe1(i),cohe2(i) coh1(i) = first cohesion coh2(i) = second cohesion cohe1(i) = strain level 1 cohe2(i) = strain level 2 NOTE: must have cohe1(i) < cohe2(i) Diagram: ________________ coh2 /| / | / | / | coh1 ______________/ cohe2 | cohe1 [GRAIN SIZE SENSITIVE CREEP] read(8,*)gscreep_mu(i),gscreep_m(i),gscreep_T(i), gscreep_strain1(i), gscreep_strain2(i),gscreep_a0(i), gscreep_a1(i), gscreep_a2(i), gscreep_updown_T [1] gscreep_mu(i) = viscosity at level of grainsize a0 NOTE: See Philippe's note on viscosities in SOPALE (there is a link to this from Bonny's web page). [2] gscreep_m(i) = power of grainsize dependency [m* in formula below] [3] gscreep_T(i) = temperature at which mu is (1/e) x that viscosity at a0 and T [T* in formula below] [4] gscreep_strain1(i) = S1* (see fig below) [5] gscreep_strain2(i) = S2* (see figure below) NOTE: Always use S1* < S2* [6] gscreep_a0(i) = parameter a0* [ a* in formula below] 0 [7] gscreep_a1(i) = a1* (see fig below) [8] gscreep_a2(i) = a2* (see fig below) [9] gscreep_updown_T [2] d m* d mu(T ) = mu* (a/a* ) exp (- T /T*) 0 [3] [6] d = dynamical [8] a1* _ _ _ _ _ _ _ _ _ ______________ /| / | / | / | [7] a2* _____________/ | | | | | S1* S2* [4] [5] d For case exp (- T /T*) gscreep_updown_T = +1 d For case exp (- T*/T ) gscreep_updown_T = -1 To suppress gs_creep in parallel, take gs_creep_mu = 0 ******** insert note about Rigid-Plastic calculation here?? **** read(8,*)nlcreep_eps(i),nlcreep_sc1(i),nlcreep_sc2(i), nlcreep_n(i),nlcreep_T(i),nlcreep_strain1(i), nlcreep_strain2(i),nlcreep_updown_T [1] nlcreep_eps(i) [2] nlcreep_sc1(i) [3] nlcreep_sc2(i) [4] nlcreep_n(i) [5] nlcreep_T(i) [6] nlcreep_strain1(i) [7] nlcreep_strain2(i) [8] nlcreep_updown_T See formula and diagram below... [4] . . n* l ( e / e* ) = (sigma/sc) exp -(T*/ T ) 0 | [5] [1] sc1=[2] if updown= -1 sc2=[3] [2] sc1* _____________ |\ | \ | \ | \ [3] sc2* _ _ _ _ _ _ | _ _\__________________ | | | | S1* S2* [6] [7] PEIERL'S LAW POWER-LAW BREAKDOWN read(8,*)pei_eps(i),pei_Tn(i),pei_sc(i),pei_T2(i), pei_sc(i),pei_T2(i),pei_strain1(i),pei_strain2(i) [1] pei_eps(i) = reference strain rate [2] pei_Tn(i) [use 0 to suppress in parallel] [3] pei_sc(i) [4] pei_T2(i) [5] pei_sc(i) [6] pei_T2(i) [7] pei_strain1(i) - not used [8] pei_strain2(i) - not used NOTE: E = epsilon [2] . . Tn*/T [4] E = E* (sigma/sc*) exp ( - T2* / T ) [1] [3] read(8,*)texpand(i),tdensity(i) texpand(i) = volume coefficient of thermal expansion tdensity(i) = reference temperature for density (in absolute degrees) if positive, density is rho(T) = rho (1 - alpha(T - Tref)) 0 [where Tref is tdensity(i)] if negative, density if rho(T) = - rho alpha T 0 [END OF INPUT FOR EACH MATERIAL'S MECHANICAL PROPERTIES] read(8,*)type_bc IF type_bc = 1 then type_bc = 1 for singularity model --------o--------- S read(8,*)vleft,vright,xsleft,xsright vleft = velocity at far left of singularity vright = velocity at far right of singularity xsleft } [see diagram] velocity changes linearly from vleft xsright } to vright across the zone xsleft xsright | | --+--------o--------+--- | S | This type of bc does not allow vertical velocities on the left or right boundaries; i.e. these are not roller boundaries. Use this type bc for models where the deformation does not reach to the ends of the Eulerian mesh. NOTE: this means there is no need to modify the Lagrangian grid that is outside the Eulerian mesh owing to vertical movements at the left and right boundaries. For cases where the left and/or right boundaries move vertically, use type_bc = 11. ------------------------------------------------- THE TYPE_BC = 3 FAMILY OF BOUNDARY CONDITIONS ------------------------------------------------- There are a number of boundary conditions based on type_bc=3, which have numbers starting with 3 (i.e. type_bc=30, type_bc=31, type_bc=312 and so on). They are variations on a theme, and have some basic properties in common, as illustrated by this diagram: +--------------------------------------+ |--> <--| |--> vxleft_up vxright_up <--| |--> <--| | | yleft1 -+- -+- yright1 | | | | yleft2 -+- -+- yright2 | | <--| |--> <--| vxleft_down vxright_down |--> <--| |--> <--| |--> <--| |--> <--| vybottom |--> <--| ^ ^ ^ ^ ^ ^ ^ ^ |--> | ^ | | | | | | | | ^ | | ^ | | | | | | | | | | ^ | +----+---------------------------+-----+ | | xvyleft xvyright NOTE: The bottom boundary y velocity is set to vybottom between xvyleft and xvyright, zero to the left of xvyleft, and zero to the right of xvyright. The variable decaylvybottom sets the length for the bottom boundary y velocity to decay to zero. In the Lagrangian grid below the Eulerian grid (i.e. where the y coordinate of the grid is less than 0), the y velocities are set to the bottom boundary y velocity. For now (as of June 4, 2002) we are setting the x velocities in this part of the Lagrangian grid to zero. This means that the mechanical and thermal materials in in this area of the Lagrangian grid are advected but not deformed. vxleft_up = x velocity at left boundary above yleft1 vxleft_down = x velocity at left boundary below yleft2 In the area between yleft1 and yleft2, the boundary condition makes a smooth transition between vxleft_up and vxleft_down. yleft1 = vertical position at left boundary, above which the x velocity is vxleft_up yleft2 = vertical position at left boundary, below which the x velocity is yleft_down rollleft = 1 means roller (free slip) boundary condition on left boundary = 0 means no vertical component to boundary condition on left boundary, i.e. pinned (no slip) vxright_up = x velocity at right boundary above yright1 vxright_down = x velocity at right boundary below yright2 In the area between yright1 and yright2, the boundary condition makes a smooth transition between vxright_up and vxright_down. yright1 = vertical position at right boundary, above which the x velocity is vxright_up yright2 = vertical position at right boundary, below which the x velocity is yright_down rollright = 1 means roller (free slip) boundary condition on right boundary = 0 means no vertical component to boundary condition on right boundary, i.e. pinned (no slip) vytop0 = 1 to set a rigid top (i.e. no vertical velocities on top boundary) = 0 vertical velocities allowed on top boundary vxbottom0 = 1 for no slip on bottom (no horizontal velocities on bottom) = 0 horizontal velocites allowed on bottom vybottom = y velocity at bottom boundary decaylvybottom = length scale in SI, m, over which vybottom (the bottom boundary y velocity) decays to zero adjacent to xleft and xright xvyleft = to the left of this point, bottom boundary y velocity is zero xvyright = to the right of this point, bottom boundary y velocity is zero -------------------------------- IF type_bc =3 then read(8,*)vxleft_up,vxleft_down,yleft1,yleft2,rollleft, vxright_up,vxright_down,yright1,yright2,rollright, vytop0,vxbottom0 vxleft_up As for type_bc=3 family vxleft_down As for type_bc=3 family yleft1 As for type_bc=3 family yleft2 As for type_bc=3 family rollleft As for type_bc=3 family vxright_up As for type_bc=3 family vxright_down As for type_bc=3 family yright1 As for type_bc=3 family yright2 As for type_bc=3 family rollright As for type_bc=3 family vytop0 As for type_bc=3 family vxbottom0 As for type_bc=3 family **** CAUTION **** type_bc = 3 has special (likely inconsistent) properties when used with roller side boundary conditions. In the tectonic calculations all nodes on the side boundaries above the base are rollers, therefore the surface nodes on the sides can move vertically. However, later when the surface is updated for tectonic advection, the vertical advection of the surface nodes at the sides is set to 0, thereby suppressing their advection (inconsistent). The advection of the surface nodes at the sides due to surface processes is calcuated independently and depends on the parameter nerode_ends. Note that all other members of the type_bc=3 family allow vertical advection of the surface nodes at the sides. type_bc=3 was the original version and should only be used in special circumstances. ---------------------------------------------- if(type_bc.eq.4)then read(8,*)vpunchx,vpunchy,xsleft,xsright endif if(type_bc.eq.5)then read(8,*)ppunchx,ppunchy,xsleft,xsright endif if(type_bc.eq.6)then read(8,*)ppunchx,xpunch,ppunchy,ypunch,xsleft,xsright endif if(type_bc.eq.7)then read(8,*)ppunchx endif if(type_bc.eq.8)then read(8,*)omegarot,domega,xsleft,xsright endif IF type_bc = 10 then read(8,*)vleft,vright,xsleft,xsright, ampnoise,wavlnoise,decaylnoise This is the same as type_bc=1 with the addition of sinusoidal noise to the x velocities between vleft and vright. vleft = velocity at far left of singularity vright = velocity at far right of singularity xsleft } [see diagram] velocity changes linearly from vleft xsright } to vright across the zone xsleft xsright | | --+--------o--------+--- | S | ampnoise = amplitude of noise as a velocity in SI, m/s wavlnoise = wavelength of sinusoidal noise in SI, m decaylnoise = length scale in SI, m, over which the noise decays to zero adjacent to xsleft and xsright Choice of noise amplitude. For the noise velocity to give a superimposed noise strain rate equal to alpha x mean velocity gradient (where mean velocity = epsilon xx tectonic strain rate), ampnoise should be set to: ampnoise = alpha x ( deltavx/ltrans) x ( wavlnoise/4 ) where deltavx = difference of velocity (m/s) across entire width of transition zone (vxright-vxleft) ltrans = width (m) of transition zone (xsright-xsleft) wavelnoise = wavelength of the noise (m) IF type_bc = 11 then type_bc = 11 for singularity model --------o--------- S read(8,*)vleft,rollleft, vright, rollright, xsleft,xsright vleft = velocity at far left of singularity rollleft 1 = allow left endpoint of surface to move vertically 0 = don't allow left endpoint of surface to move vertically vright = velocity at far right of singularity rollright 1 = allow right endpoint of surface to move vertically 0 = don't allow right endpoint of surface to move vertically xsleft } [see diagram] velocity changes linearly from vleft xsright } to vright across the zone xsleft xsright | | --+--------o--------+--- | S | WARNING! type_bc=11 allows vertical rolling of the the left and right boundaries, except at the base. No modification is made to the Lagrangian grid outside the Eulerian domain. This means there will be a mismatch between the Lagrangian and Eulerian meshes at boundaries where rolling has occurred. This is probably okay for cases where material is advected out of the Eulerian domain. It is also okay for boundaries where there is no advection into the domain. It is up to the user to make sure this is satisfied. IF type_bc = 12 then This is like type_bc=11, but allows x velocities to vary along side boundaries (as type_bc 3 does). read(8,*) vxleft_up, vxleft_down, yleft1, yleft2, rollleft, vxright_up, vxright_down, yright1, yright2, rollright, xsleft,xsright vxleft_up = x velocity at left boundary above yleft1 vxleft_down = x velocity at left boundary below yleft2 In the area between yleft1 and yleft2, the boundary condition makes a smooth transition between vxleft_up and vxleft_down. yleft1 = vertical position at left boundary, above which the x velocity is vxleft_up yleft2 = vertical position at left boundary, below which the x velocity is yleft_down rollleft 1 = allow left endpoint of surface to move vertically 0 = don't allow left endpoint of surface to move vertically vxright_up = x velocity at right boundary above yright1 vxright_down = x velocity at right boundary below yright2 In the area between yright1 and yright2, the boundary condition makes a smooth transition between vxright_up and vxright_down. yright1 = vertical position at right boundary, above which the x velocity is vxright_up yright2 = vertical position at right boundary, below which the x velocity is yright_down rollright 1 = allow right endpoint of surface to move vertically 0 = don't allow right endpoint of surface to move vertically xsleft } [see diagram] velocity changes linearly from vleft xsright } to vright across the zone xsleft xsright | | --+--------o--------+--- | S | WARNING! type_bc=12 allows vertical rolling of the the left and right boundaries, except at the base. No modification is made to the Lagrangian grid outside the Eulerian domain. This means there will be a mismatch between the Lagrangian and Eulerian meshes at boundaries where rolling has occurred. This is probably okay for cases where material is advected out of the Eulerian domain. It is also okay for boundaries where there is no advection into the domain. It is up to the user to make sure this is satisfied. IF type_bc = 30 then read(8,*)vxleft_up,vxleft_down,yleft1,yleft2,rollleft, vxright_up,vxright_down,yright1,yright2,rollright, vytop0, vxbottom0,vybottom,decaylvybottom, ampnoise,wavlnoise,decaylnoise,xsleft,xsright type_bc = 30 is like type_bc=3, but with added sinusoidal noise in the x velocity along the base. NOTE that choosing this bc means you have a no-slip bottom boundary condition, no matter what vxbottom0 is set to. vxleft_up As for type_bc=3 family vxleft_down As for type_bc=3 family yleft1 As for type_bc=3 family yleft2 As for type_bc=3 family rollleft As for type_bc=3 family vxright_up As for type_bc=3 family vxright_down As for type_bc=3 family yright1 As for type_bc=3 family yright2 As for type_bc=3 family rollright As for type_bc=3 family vytop0 As for type_bc=3 family vxbottom0 Value of vxbottom0 is ignored here, since choosing this boundary condition means you have a no-slip bottom boundary. vybottom As for type_bc=3 family decaylvybottom As for type_bc=3 family xvyleft As for type_bc=3 family xvyright As for type_bc=3 family ampnoise = amplitude of noise as a velocity in SI, m/s wavlnoise = wavelength of sinusoidal noise in SI, m decaylnoise = length scale in SI, m, over which the noise decays to zero adjacent to xsleft and xsright xsleft= left limit in SI, m, for the added basal noise function and for y velocity on bottom boundary xsright= right limit in SI, m, for the added basal noise function and for y velocity on bottom boundary if(type_bc.eq.31)then read(8,*)vxleft_up,vxleft_down,yleft1,yleft2,rollleft, & vxright_up,vxright_down,yright1,yright2,rollright, & vytop0, vxbottom0,vybottom,decaylvybottom, & ampnoise,decaylnoise,xsleft,xsright type_bc = 31 is like type_bc=3, but with added white noise in the x velocity along the base. NOTE that choosing this bc means you have a no-slip bottom boundary condition, no matter what vxbottom0 is set to. vxleft_up As for type_bc=3 family vxleft_down As for type_bc=3 family yleft1 As for type_bc=3 family yleft2 As for type_bc=3 family rollleft As for type_bc=3 family vxright_up As for type_bc=3 family vxright_down As for type_bc=3 family yright1 As for type_bc=3 family yright2 As for type_bc=3 family rollright As for type_bc=3 family vytop0 As for type_bc=3 family vxbottom0 Value of vxbottom0 is ignored here, since choosing this boundary condition means you have a no-slip bottom boundary. vybottom As for type_bc=3 family decaylvybottom As for type_bc=3 family xvyleft As for type_bc=3 family xvyright As for type_bc=3 family ampnoise = amplitude of noise as a velocity in SI, m/s decaylnoise = length scale in SI, m, over which the noise decays to zero adjacent to xsleft and xsright xsleft= left limit in SI, m, for the added basal noise function and for y velocity on bottom boundary xsright= right limit in SI, m, for the added basal noise function and for y velocity on bottom boundary IF type_bc = 312 then read(8,*)vxleft_up,vxleft_down,yleft1,yleft2,rollleft, vxright_up,vxright_down,yright1,yright2,rollright, vytop0, vxbottom0,vybottom,decaylvybottom, ampnoise,decaylnoise,xsleft,xsright type_bc = 312 is like type_bc = 3, but with added white noise in the y velocity along the base. NOTE that choosing this bc means you can have a free-slip bottom boundary condition - an improvement from type_bc.eq.31! vxleft_up As for type_bc=3 family vxleft_down As for type_bc=3 family yleft1 As for type_bc=3 family yleft2 As for type_bc=3 family rollleft As for type_bc=3 family vxright_up As for type_bc=3 family vxright_down As for type_bc=3 family yright1 As for type_bc=3 family yright2 As for type_bc=3 family rollright As for type_bc=3 family vytop0 As for type_bc=3 family vxbottom0 As for type_bc=3 family vybottom As for type_bc=3 family decaylvybottom As for type_bc=3 family xvyleft As for type_bc=3 family xvyright As for type_bc=3 family ampnoise = amplitude of noise as a velocity in SI, m/s decaylnoise = length scale in SI, m, over which the noise decays to zero adjacent to xsleft and xsright xsleft= left limit in SI, m, for the added basal noise function and for y velocity on bottom boundary xsright= right limit in SI, m, for the added basal noise function and for y velocity on bottom boundary IF type_bc = 32 then read(8,*)vxleft_up,vxleft_down,yleft1,yleft2,rollleft, vxright_up,vxright_down,yright1,yright2,rollright, vytop0, vxbottom0,vybottom,decaylvybottom, ampnoise,wavlnoise,decaylnoise,xsleft,xsright, noisenodey1,noisenodey2 type_bc = 32 is like type_bc = 3, but with added sinusoidal noise in the x body force along the nodes in rows noisenodey1 to noisenodey2 (inclusive). vxleft_up As for type_bc=3 family vxleft_down As for type_bc=3 family yleft1 As for type_bc=3 family yleft2 As for type_bc=3 family rollleft As for type_bc=3 family vxright_up As for type_bc=3 family vxright_down As for type_bc=3 family yright1 As for type_bc=3 family yright2 As for type_bc=3 family rollright As for type_bc=3 family vytop0 As for type_bc=3 family vxbottom0 As for type_bc=3 family vybottom As for type_bc=3 family decaylvybottom As for type_bc=3 family xvyleft As for type_bc=3 family xvyright As for type_bc=3 family ampnoise = amplitude of noise as a velocity in SI, m/s wavlnoise = wavelength of sinusoidal noise in SI, m decaylnoise = length scale in SI, m, over which the noise decays to zero adjacent to xsleft and xsright xsleft= left limit in SI, m, for the added noise function on rows noisenodey1 to noisenodey2 and for y velocity on bottom boundary xsright= right limit in SI, m, for the added noise function on rows noisenodey1 to noisenodey2 and for y velocity on bottom boundary noisenodey1 = top row number where noise is to be added to x body force noisenodey2 = bottom row number where noise is to be added to x body force IF type_bc = 33 then read(8,*)vxleft_up,vxleft_down,yleft1,yleft2,rollleft, vxright_up,vxright_down,yright1,yright2,rollright, vytop0, vxbottom0, vybottom,decaylvybottom, ampnoise,decaylnoise,xsleft,xsright, noisenodey1, noisenodey2 type_bc = 33 is like type_bc = 3, but with added white noise in the x body force along the nodes in rows noisenodey1 to noisenodey2 (inclusive). vxleft_up As for type_bc=3 family vxleft_down As for type_bc=3 family yleft1 As for type_bc=3 family yleft2 As for type_bc=3 family rollleft As for type_bc=3 family vxright_up As for type_bc=3 family vxright_down As for type_bc=3 family yright1 As for type_bc=3 family yright2 As for type_bc=3 family rollright As for type_bc=3 family vytop0 As for type_bc=3 family vxbottom0 As for type_bc=3 family vybottom As for type_bc=3 family decaylvybottom As for type_bc=3 family xvyleft As for type_bc=3 family xvyright As for type_bc=3 family ampnoise = amplitude of noise as a velocity in SI, m/s decaylnoise = length scale in SI, m, over which the noise decays to zero adjacent to xsleft and xsright xsleft = left limit in SI, m, for the added noise function on rows noisenodey1 to noisenodey2 and for y velocity on bottom boundary xsright = right limit in SI, m, for the added noise function on rows noisenodey1 to noisenodey2 and for y velocity on bottom boundary noisenodey1 = top row number where noise is to be added to x body force noisenodey2 = bottom row number where noise is to be added to x body force [Eulerian remeshing] read(8,*)ksplitx1,tsmooth,ismooth ksplitx1 = number of subintervales per Eulerian delta X at which markers are embedded for top surface tracking (use 2, or < 4) tsmooth = number of timesteps between two numerical smoothings of the top surface, i.e. smoothing is done every tsmooth timesteps. NOTE: in the timestep when smoothing is done, smoothing is done ncyc2 times (see subroutine vylifts). ** CAUTION *** Surface smoothing may not be consistent with all types of tectonic and surface processes. In earlier versions surface smoothing is applied to the vertically advected Eulerian tectonic surface and then surface and other processes are applied to the smoothed surface. in the nov21-02 version of the code, surface smoothing is applied to the final surface after the effects of advection owing to the surface processes have been added. Numerical surface smoothing is best avoided if possible. Instead, it is preferable to use the diffusion option in the surface processes - see subroutine erode. If smoothing is all that is required (not true surface processes) ensure that the material properties of any injected Lagrangian particles that correspond to sediment are the same as those of the material that comprises the surface that is being smoothed. ismooth = number of smoothing passes per smoothing timestep (= 1) read(8,*)vscale,verror,verror1,miterv,irescale,kinitvis, force_iter1, force_iter, force_continue, mech_trouble_iter vscale = velocity scale for scaling velocity errors in Picard iterations. Convergence is defined to be achieved when maxdv/vscale .le. verror where maxdv is the greatest change in velocity (in either x or y direction) since the previous iteration. verror = relative error for timesteps > 1 in Picard iteration verror1 = relative error for timestep 1 in Picard iterations miterv = max number of iterations in Picard iterations. If convergence is not achieved by miterv iterations, the model will stop unless force_continue is set greater than 0. irescale = scheme used for updating viscosities 0 = calculate viscosities according to material properties if stress above yield viscosities are rescaled -> viscosity rescaling created by strain rate (see earlier notes) kinitvis = max number of iterations per step before the viscosity is reinitialized and computation is started from scratch. Should always check by setting kinitvis = 1. force_iter1 = number of iterations to force for the first timestep. This number of iterations will always be performed even if convergence is achieved earlier. If convergence has not been reached after the forced iterations, the code will continue to iterate until convergence is reached or until the maximum iterations (miterv) is reached. 0 = don't force a minumum number of iterations (normal behaviour) force_iter = number of iterations to force for timesteps > 1 in the Picard iterations. This number of iterations will always be performed even if convergence is achieved earlier. If convergence has not been reached after the forced iterations, the code will continue to iterate until convergence is reached or until the maximum iterations (miterv) is reached. 0 = don't force a minimum number of iterations. (normal behaviour) force_continue = flag to indication whether to force continued timesteps when Picard divergence occurs. 0 = if Picard divergence occurs, stop the model (normal behaviour) 1 = if Picard divergence occurs, don't stop the model but continue to next timestep anyway. USE WITH CAUTION! ONLY TO BE USED IN SPECIAL CASES! mech_trouble_iter = iteration number at which to write an output file. This output file is given frame number 99, which means that its name will be something like "g01_p00_f99_o" or "g02_p00_f99_o". This file is useful for debugging problem models where convergence is not reached on the last timestep and the restart file doesn't give any clues why. If mech_trouble_iter is not specified (as in earlier versions of SOPALE1_i), it will be set to a default value of 50 iterations. read(8,*)ncyc2,color2op,lagtest ncyc2 = number of cycles in Lag/Eul remeshing/tracking per step (use 2) includes top surface recycling (0 gives no tracking) color2op = 0 -> Lagrangian particles color the Eulerian grid when computation initialized (Microfem-like process) 1 -> Eulerian grid colors Lagrangian particles initially 2 -> Philippe's code says this means "nonmixedelementcoloring grid 1 with boxes/strings" 10 -> read intial Lagrangian particle colors (materials) from file lagcolors_i lagtest - use 0 1 -> tracking with uniform velocity field, useful for benchmarks read(8,*)ucompute,tcompute,ale_compute,dt_compute,benchmark, tbottom,ttop,e_temperature,cfl, lambda_blend1, lambda_blend2 ucompute = 1 -> compute a velocity using an FE computation = 0 -> no FE computation, computation of another sort tcompute = 1 -> compute a temperature field using an FE computation = 0 -> no FE computation 0 for no thermo-mech ale_compute = 1 -> ALE computation with remeshing and computation = 0 -> no ALE computation /no remeshing dtcompute = 0 -> not going to recompute dt = 1 -> CFL recomputation of dt tbottom = bottom temperature (Absolute temperature, not Celsius) Linear geotherm top to bottom, no radioactive heating ttop = top temperature (absolute temperature) e_temperature = 0 -> Eulerian temps that control material props are elemental (one value only) = 1 -> Eulerian temps that control material props are nodal Should use 0 when temperatures are not computed, i.e. they are just advected by markers, i.e. Lag node ----------------> Eulerian element | ^ | | v | Average of lag nodes => ----| When no lag nodes in element, keep value Temperatures are not updated Entity is 'per cell' level, i.e. each element has only one property and properties are not specified at the Gauss point level cfl = used to automatically adjust timestep length (dt) when dt_compute = 1. Timestep length is reset so that the cfl takes value given in cfl. lambda_blend1, lambda_blend2 = interval during the first timestep where you want to "blend in" the value of lambda_pore_prs (fluid pore pressure ratio) in the calculation of plastic yielding (splas). In timestep 1, starting at iteration number lambda_blend1, the value of lambda_pore_prs is gradually increased from zero until it reaches its full value at iteration number lambda_blend2. This technique avoids problems with negative values of plastic yielding in the first timestep, when the model has not yet stabilized. Normally used in conjunction with plas_fluid_method (see above). NOTE: You should make sure that force_iter1 is at least as large as lambda_blend2. read(8,*)steady,analy,esteady steady = 1 for steady state computation = 0 for non-steady-state computation (steady state implies dt can be 0) analy = obsolete, set to 0 esteady = only used when steady = 1 = the error level you want to wait for until you reach steady state *** Start thermal section (only needed if tcompute = 1) **** read(8,*)maxter,dtempm,ther_advection,th_gamma_time,th_cn,tempsmooth maxter = maximum error for thermal iterations For steady=0, this is a relative error For steady=2, this is an absolute error dtempm = temperature error (units are temperature). If temperature error more than this, continue iterating. MAY NOT BE USED (value is set in code) ther_advection = 1 there is advection of the temperature field = 0 there is no advection of the temperature field th_gamma_time = numerical parameter for integration in time Choose between .5 and 1, typically 1 Trapezoidal schemes for time integration th_cn = Krank-Nicholson parameter Use 1 tempsmooth = timestep interval between successive smoothings of temperature field read(8,*)type_bct type_bct = boundary condition type for thermal problem during the model run (not initialization). See type_bct0 for initial conditions in dmtm. 1 = convection 2 = convection 3 = for subduction models for deep model thermo-mechanical (dmtm) If type_bct = 3, next line defines bct3 parameters (11 of them). This line is not needed if type_bct is not 3. bct3(1) = top temperature (degrees K) [see tshift in the code] bct3(2) = 0 bottom is imposed temperature, value bct3(3) = 1 bottom is imposed flux, value bct3(3) bct3(3) = temperature or flux at base (depending on value of bct3(2) bct3(4) = 0 On left boundary, impose initial geotherm above bct3(5) and bct3(7) for flux below bct3(5) = 1 On left boundary, impose bct3(6) as temperature above bct3(5) and bct3(7) as flux below bct3(5) = 2 On left boundary, impose bct3(6) as temperature above bct3(5) and initial geotherm below bct3(5) bct3(5) = z Position on left boundary about which the boundary condition changes (see bct3(4)) - distance in metres or user units measured from origin upwards bct3(6) = T Temperature on left boundary when bct3(4) =1 or 2 Not used when bct3(4) = 0. bct3(7) = flux Lateral flux imposed at left boundary below bct3(5) when bct3(4) = 0 or 1. Not used when bct3(4) = 2. bct3(8) = 0 On right boundary, impose initial geotherm above bct3(9) and bct3(11) for flux below bct3(9) = 1 On right boundary, impose bct3(10) as temperature above bct3(9) and bct3(11) as flux below bct3(9) = 2 On right boundary, impose bct3(10) as temperature above bct3(9) and initial geotherm below bct3(9) bct3(9) = z Position on right boundary about which the boundary condition changes (see bct3(8)) - distance in metres or user units measured from origin upwards bct3(10) = T Temperature on right boundary when bct3(8) =1 or 2 Not used when bct3(8) = 0. bct3(11) = flux Lateral flux imposed at right boundary below bct3(9) when bct3(8) = 0 or 1. Not used when bct3(8) = 2. read(8,*)th_dt,cfl,vzgiven th_dt = timestep for thermal code Is computed internally when steady=0 and dtcompute=1 (steady=0 means non steady state calculation) cfl NOTE: cfl is re-read; the value used is the last value taken vzgiven = only used for benchmarks 20, 21, 22, 23 not used in deep thermo-mechanical models read(8,*)ncolor2t ncolor2t = number of thermal materials For each thermal material (1 to ncolor2t), have the following 5 lines: read(8,*)thermod1_conductivity(i) = K, thermal conductivity units W/(m K) or user units read(8,*)thermod1_rau0(i) = Thermal density Density used by thermal code units kg/m**3 or user units read(8,*)thermod1_cp(i) = Specific heat units J/(kg K) or user units read(8,*)thermod1_fheating(i) = frictional heating coefficient Not active, use 0 read(8,*)thermod1_power_permass(i) = internal heat production (e.g. radiactive heat generation) units watts/kg (NOTE: need to convert from A in units of watts/m**3 to H in units of watts/kg. A/rau0 = H) *** End thermal section (only needed if tcompute = 1) **** READ(8,*) OUTFIL OUTFIL = name of output file in which blkfct will write diagnostic info Per linear SPD system (i.e. per blkfct system) need 4 lines for each of MECH and THERMAL SYSTEM if they are needed READ(8,*) MATFIL MATFIL = file name containing profile of the matrix to be solved Can have 0, 1, or 2 systems to solve These are all in "matrices subdirectory" and are generated by program pmain.f in subdirectory "sources". READ(8,*) ICASE ICASE = 1 natural ordering = 2 mmd ordering READ(8,*) CACHSZ CACHSZ = cachesize in kilobytes Try 64, can be flexible READ(8,*) LEVEL LEVEL = loop unrolling level Try 8 *** Begin thermal section (only needed if tcompute = 1) **** READ(8,*) MATFIL (for thermal systems) MATFIL = file name containing profile of the matrix to be solved Can have 0, 1, or 2 systems to solve These are all in "matrices subdirectory" and are generated by program pmain.f in subdirectory "sources". READ(8,*) ICASE (for thermal systems) ICASE = 1 natural ordering = 2 mmd ordering READ(8,*) CACHSZ (for thermal systems) CACHSZ = cachesize in kilobytes Try 64, can be flexible READ(8,*) LEVEL (for thermal systems) LEVEL = loop unrolling level Try 8 *** End thermal section (only needed if tcompute = 1) *** read(8,*)ibinread ibinread = 1 binary read/writes to outputs = 0 formatted read/writes (may not work) NOTE: the following filenames are compulsory in terms of number of characters, i.e. the filenames cannot have any more or less characters than the format specified in the read statement. Also, the first 13 characters of the filename are assumed to specify the directory (including the trailing "/"). read(8,'(a26)')filename1 filename1 = name of restart file for Eulerian grid read(8,'(a26)')filename2 filename2 = name of restart file for Lagrangian grid read(8,'(a33)')filenameg1 filenameg1 = name of file to save timestep information read(8,'(a26)')filenameg3 filenameg3 = name of file containing grid type and topology used by mozart_plot.x for Eulerian read(8,'(a26)')filenameg4 filenameg4 = name of file containing grid type and topology used by mozart_plot.x for Lagrangian *** Mechanical material (colour) boxes section **** NOTE: THIS INPUT DEPENDS ON color2op. If color2op = 0, we will read these lines read(8,*)ncolor2 ncolor2 = number of Lagrangian mechanical material boxes defined on lag grid Then for each box, read the 4 vertices (in metres or user units) of the box in counterclockwise order, i.e.: read(8,*)xx1,xx2,xx3,xx4 xx1 = x coordinate of top left corner of box xx2 = z coordinate of top left corner of box xx3 = x coordinate of bottom left corner xx4 = z coordinate of bottom left corner read(8,*)xx5,xx6,xx7,xx8 xx5 = x coordinate of bottom right corner xx6 = z coordinate of bottom right corner xx7 = x coordinate of top right corner xx8 = z coordinate of top right corner read(8,*)newcol2 newcol2 = corresponding material number for material particles in this box If color2op = 1 then define the layers on Eulerian grid : read(8,*)nelayers nelayers = number of layers For each layer define the last row of the Eulerian grid for this layer, counting top to bottom read(8,*)(layerse(i),i=1,nelayers) layerse(i) are the last row on the Eulerian grid for each of the nelayers Each layer is automatically given the mechanical material (colour) number in order of the layers, i.e. layer 1 is given material number 1, layer 2 is material 2, etc. Round off error determines box. NOTE: Slight perturbation on lag grid = 0.333 x 10e-6 in x Mixing model determines material property number to be taken by an Eulerian element during regridding when there are Lagrangian particles with different material property numbers in the Eulerian element. The mixing model either uses the minimum material property number (color number), or the material property number of the last particle found (*** needs to be checked ***). *** Thermal material (colour) boxes section *** DMTM version note: You need to define at least one thermal material box, even if you're not doing thermal calculations. If you are doing a mechanical-only model, define a dummy thermal material box which covers your entire domain and consists of thermal material 1. (You don't have to define a thermal material). read(8,*)ncolor2 ncolor2 = number of Lagrangian thermal material boxes defined on lag grid Then for each box, read the 4 vertices (in metres or user units) of the box in counterclockwise order, i.e.: read(8,*)xx1,xx2,xx3,xx4 xx1 = x coordinate of top left corner of box xx2 = z coordinate of top left corner of box xx3 = x coordinate of bottom left corner xx4 = z coordinate of bottom left corner read(8,*)xx5,xx6,xx7,xx8 xx5 = x coordinate of bottom right corner xx6 = z coordinate of bottom right corner xx7 = x coordinate of top right corner xx8 = z coordinate of top right corner read(8,*)newcol2 newcol2 = corresponding thermal material number for material particles in this box *** End of thermal material (colour) boxes section *** read(8,*)dmtm dmtm = 1 flag for deep thermo-mechanical model dmtm = 0 if you are doing mechanical model only Only if dmtm = 1 and rres is NOT 1, then read the following initial thermal boundary conditions: read(8,*)type_bct0 type_bct0 = initial thermal boundary condition type = 3 (only value currently implemented) for deep model thermo-mechanical (dmtm) NOTE: type_bct0 MUST BE 3!! The code assumes that values of bct30(i) [see following] are applied before trying to solve the initial thermal matrix. (We tried type_bct0=0 and it there was a "MATRIX IS NOT POSITIVE DEFINITE" error.) If type_bct0 = 3, next line defines bct30 parameters (11 of them). read(8,*)(bct30(i),i=1,11) bct30(1) = top temperature (degrees K) [see tshift in the code] bct30(2) = 0 bottom is imposed temperature, value bct30(3) = 1 bottom is imposed flux, value bct30(3) bct30(3) = temperature or flux at base (depending on value of bct30(2) bct30(4) = 0 On left boundary, impose initial geotherm above bct30(5) and bct30(7) for flux below bct30(5) = 1 On left boundary, impose bct30(6) as temperature above bct30(5) and bct30(7) as flux below bct30(5) = 2 On left boundary, impose bct30(6) as temperature above bct30(5) and initial geotherm below bct30(5) bct30(5) = z Position on left boundary about which the boundary condition changes (see bct30(4)) - distance in metres or user units measured from origin upwards bct30(6) = T Temperature on left boundary when bct30(4) =1 or 2 Not used when bct30(4) = 0. bct30(7) = flux Lateral flux imposed at left boundary below bct30(5) when bct30(4) = 0 or 1. Not used when bct30(4) = 2. bct30(8) = 0 On right boundary, impose initial geotherm above bct30(9) and bct30(11) for flux below bct30(9) = 1 On right boundary, impose bct3(10) as temperature above bct30(9) and bct30(11) as flux below bct30(9) = 2 On right boundary, impose bct30(10) as temperature above bct30(9) and initial geotherm below bct30(9) bct30(9) = z Position on right boundary about which the boundary condition changes (see bct30(8)) - distance in metres or user units measured from origin upwards bct30(10) = T Temperature on right boundary when bct30(8) =1 or 2 Not used when bct30(8) = 0. bct30(11) = flux Lateral flux imposed at right boundary below bct30(9) when bct30(8) = 0 or 1. Not used when bct30(8) = 2. NOTE no geotherm Full calculation of diffusion/production thermal equilibrium Note need for consistent temperature input - use Kelvin read(8,*)ntperts ntperts = number of boxes for temperature perturbation For each box, read in coordinates of corners and temperature perturbation value (tpert) do k=1,ntperts read(8,*)xx1,xx2,xx3,xx4 read(8,*)xx5,xx6,xx7,xx8 read(8,*)tpert =========================================================================== File: sopale1_plus_i This file contains parameters to define the Eulerian grid, the properties of the temperature field on the Eulerian grid, the Lagrangian grid, and the properties of the temperature field on the Lagrangian grid. ** NOTE ** (0,0) is defined to be the origin of the Eulerian grid. This means that Eulerian nodes will have only positive coordinates. The Lagrangian grid, however, may have negative coordinates. First section: EULERIAN grid read(78,*)xor,lhoriz,nlayers,nbzeremesh xor = x origin, in user units [MUST BE ZERO] lhoriz = grid width, in user units nlayers = number of layers with different spacing in grid nbzeremesh = row number in Eulerian grid, above which grid is uniformly dilated (stretched) during Eulerian regridding to account for changes in the top or bottom surfaces. Generally chosen to match a change in the vertical refinement of the Eulerian grid. See subroutine zeremesh. When the whole grid is vertically dilated nbzeremesh = ny (total number of nodes vertically in the Eulerian grid.) read(78,*)(yl(i),i=1,nlayers+1) yl(i) = depth (in user units) of the layer interfaces (i.e. depths of the boundaries between the layers). These are defined from top to bottom, where the bottom is at depth 0. read(78,*)(nyl(i),i=1,nlayers) nyl(i) = number of elements in layer i Note that the layers are numbered from top to bottom. Next section: EULERIAN THERMAL grid read(78,*) nlt nlt = number of layers in the thermal grid read(78,*) ytop(i),ylthick(i),rada(i),cond(i),temptopK(i),qbase(i) ytop = position of top interface of layer in model coordinate system ylthick = thickness of layer rada = radioactive heat generation, A, ie production per unit volume cond = thermal conductivity of layer temptopK = temperature of upper interface of the layer (Kelvin) qbase = vertical heat flux into the base of the layer Next section: LAGRANGIAN grid read(78,*)xor,lhoriz,nlayers xor = x origin, in user units lhoriz = grid width, in user units nlayers = number of layers with different spacing in grid read(78,*)(yl(i),i=1,nlayers+1) yl(i) = depth (in user units) of the layer interfaces (i.e. depths of the boundaries between the layers). These are defined from top to bottom, where the bottom is at depth 0. read(78,*)(nyl(i),i=1,nlayers) nyl(i) = number of elements in layer i Note that the layers are numbered from top to bottom. Next section: LAGRANGIAN THERMAL grid read(78,*) nlt nlt = number of layers in the thermal grid read(78,*) ytop(i),ylthick(i),rada(i),cond(i),temptopK(i),qbase(i) ytop = position of top interface of layer in model coordinate system ylthick = thickness of layer rada = radioactive heat generation, A, ie production per unit volume cond = thermal conductivity of layer temptopK = temperature of upper interface of the layer (Kelvin) qbase = vertical heat flux into the base of the layer =========================================================================== File: erode_i read(76,*)erosion_model erosion_model = type of erosion model you want to use = 0 for no erosion = 1 for combinations of total erosion and total deposition depending on nerode_type = 2 for slope-dependent erosion/deposition = 3 for diffusive erosion (added by Susanne Buiter) read(76,*)dt_e dt_e = timestep length in seconds (must be equal to dt in SOPALE1_i) NOTE: this is now obsolete. The SOPALE code calculates the appropriate dt_e (it uses dt_adv) instead of relying on this value. However, you must still put dt_e in the erode_i input file. read(76,*)erode_p(1) If erosion_model = 2 (slope dependent erosion) then erode_p(1) = erosion rate per second when the surface slope = 1 If erosion_model = 3 (diffusive erosion) then erode_p(1) = diffusion coefficient k (in m^2/s) read(76,*)nerode_type, nerode_ends nerode_type = 1 erosion only = 2 erosion and deposition (may be unstable, gives element-by-element oscillations in Eulerian grid) = 3 deposition only nerode_ends = 0 No erosion of surface endpoints This is the recommended value when deformation in the model will not propagate to the boundaries = 1 Erodes points 1 and nx by same amount as points 2 and nx-1, respectively. This gives no-slope boundaries at the edges for surface processes. = 2 Erodes points 1 and nx using local slope NOT CODED YET read(76,*)er_stablax, method_slope_avg er_stablax = Lax stabilization factor The larger this value is, the more stabilization is applied. method_slope_avg = 1 calculate average slope by taking the absolute values of the slopes and then averaging = 2 calculate average slope by taking the slope between points 1 and 3. If this line (containing er_stablax and method_slope_avg) is missing from erode_i, the code should set er_stablax to 1 and method_slope_avg to 1. Look in the first few lines of the standard output (so.xout) to see what the code is using for er_stablax and method_slope_avg. =========================================================================== File: initsurf_i This file contains the parameters to define an initial surface topography. The surface in initsurf_i is used to modify the grids as defined in SOPALE1_i and sopale1_plus_i. If initsurf_i is not found by SOPALE, it will assume that the grid as defined by SOPALE1_i and sopale1_plus_i does not need any further modifications. NOTE: There is currently (as of July 16, 2002) a bug in SOPALE which causes the model to die at the beginning if we define a surface which rises above the rectangular grid as defined in SOPALE1_i. Until we fix this bug, you should make sure your initial surface topology lies completely below the height of the Eulerian grid. read(10,*)refh refh = reference height for erosion (if using erosion, the surface above this height will be eroded, surface below this height will be filled in with deposition) Next section: EULERIAN grid read(10,*)nbinitmesh nbinitmesh = row number in grid (counting from the top), above which the vertical grid spacing is dilated or compressed to accomodate the surface topology; at this row and below, the grid spacing remains whatever was defined in sopale1_plus_i. NOTE: make sure you choose nbinitmesh to be further down than the lowest point of your defined surface topography. read(10,*)ntopdat ntopdat > 0 = number of points in surface function (points to be defined in next ntopdat rows) which are joined with straight-line segments to define the surface topography = -1 define a sinuoisal surface topography Next lines: if ntopdat > 0 then the next ntopdat lines are: read(10,*)topdat(i,1),topdat(i,2) topdat(i,1) = x coordinate of point i topdat(i,2) = y coordinate of point i if ntopdat = -1, then read(10,*) xsleft_surf, xsright_surf, amp_surf, wavl_surf, decayl_surf,ymid_sin xsleft_surf - left endpoint of sinuosoidal surface topography xsright_surf - right endpoint of sinuosoidal surface topography amp_surf - amplitude of sine wave wavl_surf - wavelength of sine wave decayl_surf - length scale over which the amplitude of the sine wave decays to zero adjacent to xsleft_surf and xsright_surf ymid_sin - y position of midline of sine wave. Choose this value so that ymid_sin + amp_surf is equal to or less than the height of the Eulerian grid. Next section: LAGRANGIAN grid (same as above, but applied to the Lagrangian grid) NOTE: Be careful that your topology is defined for the entire length of the grid. Otherwise you may get a "cliff" where the surface rises abruptly to the original height of the grid. read(10,*)nbinitmesh nbinitmesh = row number in grid (counting from the top), above which the vertical grid spacing is dilated or compressed to accommodate the surface topology; at this row and below, the grid spacing remains whatever was defined in sopale1_plus_i. NOTE: make sure you choose nbinitmesh to be further down than the lowest point of your defined surface topography. read(10,*)ntopdat ntopdat > 0 = number of points in surface function (points to be defined in next ntopdat rows) which are joined with straight-line segments to define the surface topography = -1 define a sinuoisal surface topography Next lines: if ntopdat > 0 then the next ntopdat lines are: read(10,*)topdat(i,1),topdat(i,2) topdat(i,1) = x coordinate of point i topdat(i,2) = y coordinate of point i if ntopdat = -1, then read(10,*) xsleft_surf, xsright_surf, amp_surf, wavl_surf, decayl_surf,ymid_sin xsleft_surf - left endpoint of sinuosoidal surface topography xsright_surf - right endpoint of sinuosoidal surface topography amp_surf - amplitude of sine wave wavl_surf - wavelength of sine wave decayl_surf - length scale over which the amplitude of the sine wave decays to zero adjacent to xsleft_surf and xsright_surf ymid_sin - y position of midline of sine wave. Choose this value so that ymid_sin + amp_surf is equal to or less than the height of the Eulerian grid.