read_input_file.f90

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              READ_INPUT_FILE    Feb. 2007                                    |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

subroutine read_input_file (params,threadinfo,material0,mat,             &
                            surface,boxes,sections,cube_faces)         

!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|

! subroutine where all default values for input parameters are defined
! it is also where the input file 'input.txt' is read and the corresponding
! parameters are set

! nstep is number of time steps
! ns is number of surface
! dt is time step length
!  (note that if dt is negative, it will be replaced by a dynamically
!   determined value derived from Courant's condition and Courant's
!   parameter called "courant")
! nmat is number of material
! material0 is material number for stuff above 1st surface
!   or reference material when there is no surface
! mat:material properties, including, density,viscosity and compressibility
!     penalty parameters; ms is dimensioned (0:nmat)
!     material 0 is void
! leveluniform_oct is level of uniform discretization (power of 2)
! levelmax_oct is maximum discretization level
! levelcut is maximum level (within a leave) used to estimate integrals
!    in divfem approximation
! levelapprox is maximum level to use to estimate positive volume
!    beyond levelcut
! noctreemax is maximum size of octrees when they are created
! penalty is penalty parameters for linear constraints arising from
!    bad faces in octree discretization (should be large)
! tempscale is temperature scaling parameter
! strain_ratio is used for refinement (if=1 no refinement;
!    if=0 uniform refinement) down to levelmax_oct
! istrain_refine: two algortihms are allowed for the refinement based on the strain ratio
!    one is based on the maximum difference between any two components of velocity
!    inside an element (istrain_refine=0) the other is based on the norm of the
!    velocity gradient (istrain_refine=1). Default is istrain_refine=1
! courant is ratio of courant conditoin used for moving particles (<1)
! stretch is the maximum allowed increase in linear length between two initially
!    adjacent particles on any surface; when this stretch is achieved, a new
!    particle is inserted on the surface, half-way along the stretched edge
!    stretch is specific to each surface
! anglemax is the maximum allowed angle between two adjacent normals; when the
!    angle is reached, a new point is injected
!    anglemax is specific to each surface
!    default value is 1.d0
! surface are surface structure
!    for each surface, one needs to define a levelt, itype, material and fnme.
!    levelt is the inital level for the particles on the surface; to be accurate
!    and avoid wholes in the surface during definition of the lsf, one should use
!    levelt=levelmax_oct+1 for all surfaces as a minimum value; itype should be 1
!    for foldable surfaces or 0 for nonfoldable surfaces; material is the material
!    type refering to the table of material available (max nmat); fnme is the name
!    of the file containing the geometry of the particles defining the surface
! npmin is the minimum number of particles in the strain cloud per element
! npmax is the maximum number of particles in the strain cloud per element
!       at levelmax_oct level (smallest possible elements)
! nonlineariter is number of iterations in nonlinear analysis
!   if nonlineariter is 0 no nonlinear iterations are performed (linear analysis)
!   if nonlineariter is positive, nonlineariter is the number of iterations performed
!                         regardless of convergence of solution
!   if nonlineariter is negative, -nonlineariter is maximum number of iterations allowed
!                         to reach convergence as determined by tol
! tol is relative tolerance (duvw/uvw) to achieve convergence
! criterion is criterion used to define the octree in the vicinity of the
!    sufaces; criterion 1 corresponds to imposing that all leaves that contain at
!    least one particle of any surface is at levelmax_oct; criterion 2 corresponds
!    to imposing that discretization is proportional to the curvature of the
!    surface; curvature is calculated from the local divergence of the normals.
!    the criterion is specific to each surface (default is 2)
! anglemaxoctree is only defined for criterion 2; it is the maximum allowable angle
!    between two adjacent normals; if the angle is greater than anglemaxoctree, the local
!    octree leaves are forced to be at level levelmax_oct; otherwise they are
!    proportionally larger (smaller levels) (default is 10)
!    anglemaxoctree is specific to each surface
!   (be aware that these files are enormous...)
! niter_move is number of iterations used to move particles in an
!  implicit, mid-point algorithm
! restart is a restart flag; if irestart is not 0, the run will restart from
!    an output file given by restartfile and at step irestart+1
! ismooth is a flag to impose an additional level of smoothing after refinement
!    for the surfaces and strain rate. It ensures that no leaf is flanked by
!    other leaves diffeing by more than 1 level of refinement
!    If ismooth is 0, no smoothing is performed; if ismooth is set to 1, smoothing
!    is performed (default is 1)
! boxes is the set of nboxes box structures defined by the user where a set level
!    of discretization is imposed

!------------------------------------------------------------------------------|
!((((((((((((((((  declaration of the subroutine arguments  ))))))))))))))))))))
!------------------------------------------------------------------------------|

use threads
use definitions
use constants 

implicit none

type (parameters) params
type (thread) threadinfo
integer material0
type (material) mat(0:params%nmat)
type (sheet) surface(params%ns)
type (box) boxes(params%nboxes)
type (cross_section) sections(params%nsections)
type (face),dimension(6) :: cube_faces

!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|

integer ires,i,il,ierr,iproc,nproc
character cm*3
character(len=1) :: answer
integer iran,doru
integer flag

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

INCLUDE 'mpif.h'

call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)
!==============================================================================
!==============================================================================

params%nstep=2
if (iproc==0) call scanfile (params%infile,'nstep',params%nstep,ires)
call mpi_bcast(params%nstep,1,mpi_integer,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'nstep ',params%nstep

!=====[material properties]====================================================

material0=0
if (iproc==0) call scanfile (params%infile,'material0',material0,ires)
call mpi_bcast(material0,1,mpi_integer,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'material0 ',material0

do i=0,params%nmat
   write(cm,'(i3)') i
   il=1
   if (i.lt.100) il=2
   if (i.lt.10) il=3

   mat(i)%density=1.d0
   if (i.eq.0) mat(i)%density=0.d0
   if (iproc==0) call scanfile (params%infile,'density'//cm(il:3),mat(i)%density,ires)
   call mpi_bcast(mat(i)%density,1,mpi_double_precision,0,mpi_comm_world,ierr)
   write(threadinfo%Logunit,*) 'density mat ',i,mat(i)%density
 
   mat(i)%viscosity=1.d0
   if (i.eq.0) mat(i)%viscosity=1.d-8
   if (iproc==0) call scanfile (params%infile,'viscosity'//cm(il:3),mat(i)%viscosity,ires)
   call mpi_bcast(mat(i)%viscosity,1,mpi_double_precision,0,mpi_comm_world,ierr)
   write(threadinfo%Logunit,*) 'viscosity mat ',i,mat(i)%viscosity

   mat(i)%penalty=1.d8
   if (iproc==0) call scanfile (params%infile,'penalty'//cm(il:3),mat(i)%penalty,ires)
   call mpi_bcast(mat(i)%penalty,1,mpi_double_precision,0,mpi_comm_world,ierr)
   write(threadinfo%Logunit,*) 'penalty mat ',i,mat(i)%penalty

   mat(i)%expon=1.d0
   if (iproc==0) call scanfile (params%infile,'expon'//cm(il:3),mat(i)%expon,ires)
   call mpi_bcast(mat(i)%expon,1,mpi_double_precision,0,mpi_comm_world,ierr)
   write(threadinfo%Logunit,*) 'expon mat ',i,mat(i)%expon
   
   mat(i)%activationenergy=0.d0
   if (iproc==0) call scanfile (params%infile,'activationenergy'//cm(il:3),mat(i)%activationenergy,ires)
   call mpi_bcast(mat(i)%activationenergy,1,mpi_double_precision,0,mpi_comm_world,ierr)
   write(threadinfo%Logunit,*) 'activationenergy mat ',i,mat(i)%activationenergy
   
   mat(i)%expansion=0.d0
   if (iproc==0) call scanfile (params%infile,'expansion'//cm(il:3),mat(i)%expansion,ires)
   call mpi_bcast(mat(i)%expansion,1,mpi_double_precision,0,mpi_comm_world,ierr)
   write(threadinfo%Logunit,*) 'expansion mat ',i,mat(i)%expansion

   mat(i)%diffusivity=1.d0
   if (iproc==0) call scanfile (params%infile,'diffusivity'//cm(il:3),mat(i)%diffusivity,ires)
   call mpi_bcast(mat(i)%diffusivity,1,mpi_double_precision,0,mpi_comm_world,ierr)
   write(threadinfo%Logunit,*) 'diffusivity mat ',i,mat(i)%diffusivity

   mat(i)%heat=0.d0
   if (iproc==0) call scanfile (params%infile,'heat'//cm(il:3),mat(i)%heat,ires)
   call mpi_bcast(mat(i)%heat,1,mpi_double_precision,0,mpi_comm_world,ierr)
   write(threadinfo%Logunit,*) 'heat mat ',i,mat(i)%heat

   mat(i)%plasticity_type='No'
   if (iproc==0) call scanfile (params%infile,'plasticity_type'//cm(il:3),mat(i)%plasticity_type,ires)
   call mpi_bcast(mat(i)%plasticity_type,8,mpi_character,0,mpi_comm_world,ierr)
   write(threadinfo%Logunit,*) 'plasticity_type mat ',i,mat(i)%plasticity_type

   mat(i)%plasticity_parameters=0.d0
   mat(i)%plasticity_parameters(5)=-1.d0
   if ( trim(mat(i)%plasticity_type).ne.'No') then
       if (iproc==0) call scanfile (params%infile,'plasticity_1st_param'//cm(il:3),mat(i)%plasticity_parameters(1),ires)
       if (iproc==0) call scanfile (params%infile,'plasticity_2nd_param'//cm(il:3),mat(i)%plasticity_parameters(2),ires)
       if (iproc==0) call scanfile (params%infile,'plasticity_3rd_param'//cm(il:3),mat(i)%plasticity_parameters(3),ires)
       if (iproc==0) call scanfile (params%infile,'plasticity_4th_param'//cm(il:3),mat(i)%plasticity_parameters(4),ires)
       if (iproc==0) call scanfile (params%infile,'plasticity_5th_param'//cm(il:3),mat(i)%plasticity_parameters(5),ires)
       if (iproc==0) call scanfile (params%infile,'plasticity_6th_param'//cm(il:3),mat(i)%plasticity_parameters(6),ires)
       if (iproc==0) call scanfile (params%infile,'plasticity_7th_param'//cm(il:3),mat(i)%plasticity_parameters(7),ires)
       if (iproc==0) call scanfile (params%infile,'plasticity_8th_param'//cm(il:3),mat(i)%plasticity_parameters(8),ires)
       if (iproc==0) call scanfile (params%infile,'plasticity_9th_param'//cm(il:3),mat(i)%plasticity_parameters(9),ires)
       call mpi_bcast(mat(i)%plasticity_parameters,9,mpi_double_precision,0,mpi_comm_world,ierr)
       write(threadinfo%Logunit,*) 'plasticity_parameters mat ',i,mat(i)%plasticity_parameters
   endif

enddo

params%viscositymin=-1.d0
if (iproc==0) call scanfile (params%infile,'viscositymin',params%viscositymin,ires)
call mpi_bcast(params%viscositymin,1,mpi_double_precision,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'viscositymin ',params%viscositymin

params%viscositymax=-1.d0
if (iproc==0) call scanfile (params%infile,'viscositymax',params%viscositymax,ires)
call mpi_bcast(params%viscositymax,1,mpi_double_precision,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'viscositymax ',params%viscositymax

params%leveluniform_oct=3
if (iproc==0) call scanfile (params%infile,'leveluniform_oct',params%leveluniform_oct,ires)
call mpi_bcast(params%leveluniform_oct,1,mpi_integer,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'leveluniform_oct',params%leveluniform_oct


params%levelmax_oct=4
if (iproc==0) call scanfile (params%infile,'levelmax_oct',params%levelmax_oct,ires)
call mpi_bcast(params%levelmax_oct,1,mpi_integer,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'levelmax_oct',params%levelmax_oct


params%levelcut=2
if (iproc==0) call scanfile (params%infile,'levelcut',params%levelcut,ires)
call mpi_bcast(params%levelcut,1,mpi_integer,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'levelcut',params%levelcut


params%levelapprox=3
if (iproc.eq.0) call scanfile (params%infile,'levelapprox',params%levelapprox,ires)
call mpi_bcast(params%levelapprox,1,mpi_integer,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'levelapprox',params%levelapprox

params%calculate_temp=.true.
if (iproc==0) then
   call scanfile (params%infile,'calculate_temp',answer,ires)
   params%calculate_temp=(trim(answer)=='T')
end if
call mpi_bcast(params%calculate_temp,1,mpi_logical,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'calculate_temp',params%calculate_temp

params%ztemp=1.d0
if (iproc==0) call scanfile (params%infile,'ztemp',params%ztemp,ires)
call mpi_bcast(params%ztemp,1,mpi_double_precision,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'ztemp',params%ztemp
 
params%debug=0
if (iproc==0) call scanfile (params%infile,'debug',params%debug,ires)
call mpi_bcast(params%debug,1,mpi_integer,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'debug',params%debug

params%smoothing_type=0
if (iproc==0) call scanfile (params%infile,'smoothing_type',params%smoothing_type,ires)
call mpi_bcast(params%smoothing_type,1,mpi_integer,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'smoothing_type',params%smoothing_type

params%normaladvect=.false.
if (iproc==0) then
   call scanfile (params%infile,'normaladvect',answer,ires)
   params%normaladvect=(trim(answer)=='T')
end if
call mpi_bcast(params%normaladvect,1,mpi_logical,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'normaladvect',params%normaladvect


params%excl_vol=.false.
if (iproc==0) then
   call scanfile (params%infile,'excl_vol',answer,ires)
   params%excl_vol=(trim(answer)=='T')
end if
call mpi_bcast(params%excl_vol,1,mpi_logical,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'excl_vol',params%excl_vol

params%adaptive_tol=.false.
if (iproc==0) then
   call scanfile (params%infile,'adaptive_tol',answer,ires)
   params%adaptive_tol=(trim(answer)=='T')
end if
call mpi_bcast(params%adaptive_tol,1,mpi_logical,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'adaptive_tol',params%adaptive_tol

!=====[surface properties]=====================================================
        
do i=1,params%ns
   write(cm,'(i3)') i
   il=1
   if (i.lt.100) il=2
   if (il.lt.10) il=3
   
   surface(i)%itype=1
   if (iproc==0) call scanfile (params%infile,'itype'//cm(il:3),surface(i)%itype,ires)
   call mpi_bcast(surface(i)%itype,1,mpi_integer,0,mpi_comm_world,ierr)
   write(threadinfo%Logunit,*) 'surface ',i,'itype',surface(i)%itype
   
   surface(i)%material=1
   if (iproc==0) call scanfile (params%infile,'material'//cm(il:3),surface(i)%material,ires)
   call mpi_bcast(surface(i)%material,1,mpi_integer,0,mpi_comm_world,ierr)
   write(threadinfo%Logunit,*) 'surface ',i,'material',surface(i)%material

   surface(i)%rand=.false.
   if (iproc==0) then
      call scanfile (params%infile,'rand'//cm(il:3),answer,ires)
      surface(i)%rand=(trim(answer)=='T')
   end if
   call mpi_bcast(surface(i)%rand,1,mpi_logical,0,mpi_comm_world,ierr)
   write(threadinfo%Logunit,*) 'surface ',i,'rand',surface(i)%rand
   
   surface(i)%levelt=params%levelmax_oct+1
   if (iproc==0) call scanfile (params%infile,'levelt'//cm(il:3),surface(i)%levelt,ires)
   call mpi_bcast(surface(i)%levelt,1,mpi_integer,0,mpi_comm_world,ierr)
   write(threadinfo%Logunit,*) 'surface ',i,'levelt',surface(i)%levelt
   
   surface(i)%stretch=1.5d0
   if (iproc==0) call scanfile (params%infile,'stretch'//cm(il:3),surface(i)%stretch,ires)
   call mpi_bcast(surface(i)%stretch,1,mpi_double_precision,0,mpi_comm_world,ierr)
   write(threadinfo%Logunit,*) 'surface ',i,'stretch',surface(i)%stretch

   surface(i)%criterion=1
   if (iproc==0) call scanfile (params%infile,'criterion'//cm(il:3),surface(i)%criterion,ires)
   call mpi_bcast(surface(i)%criterion,1,mpi_integer,0,mpi_comm_world,ierr)
   write(threadinfo%Logunit,*) 'surface ',i,'criterion',surface(i)%criterion

   surface(i)%anglemax=1.d0
   if (iproc==0) call scanfile (params%infile,'anglemax'//cm(il:3),surface(i)%anglemax,ires)
   surface(i)%anglemax=surface(i)%anglemax*pi/180.d0
   call mpi_bcast(surface(i)%anglemax,1,mpi_double_precision,0,mpi_comm_world,ierr)
   write(threadinfo%Logunit,*) 'surface ',i,'anglemax',surface(i)%anglemax

   surface(i)%anglemaxoctree=1.d0
   if (iproc==0) call scanfile (params%infile,'anglemaxoctree'//cm(il:3),surface(i)%anglemaxoctree,ires)
   surface(i)%anglemaxoctree=surface(i)%anglemaxoctree*pi/180.d0
   call mpi_bcast(surface(i)%anglemaxoctree,1,mpi_double_precision,0,mpi_comm_world,ierr)
   write(threadinfo%Logunit,*) 'surface ',i,'anglemaxoctree',surface(i)%anglemaxoctree

   surface(i)%spread_surface_points=0
   if (iproc==0) call scanfile (params%infile,'spread_surface_points'//cm(il:3),surface(i)%spread_surface_points,ires)
   call mpi_bcast(surface(i)%spread_surface_points,1,mpi_integer,0,mpi_comm_world,ierr)
   write(threadinfo%Logunit,*) 'surface ',i,'spread_surface_points',surface(i)%spread_surface_points

   if (iproc==0) call scanfile (params%infile,'surface_type_'//cm(il:3),surface(i)%surface_type,ires)
   call mpi_bcast(surface(i)%surface_type,1,mpi_integer,0,mpi_comm_world,ierr)
   write(threadinfo%Logunit,*) 'surface ',i,'type',surface(i)%surface_type
 
   if (iproc==0) call scanfile (params%infile,'surface_param_01_'//cm(il:3),surface(i)%sp01,ires)
   if (iproc==0) call scanfile (params%infile,'surface_param_02_'//cm(il:3),surface(i)%sp02,ires)
   if (iproc==0) call scanfile (params%infile,'surface_param_03_'//cm(il:3),surface(i)%sp03,ires)
   if (iproc==0) call scanfile (params%infile,'surface_param_04_'//cm(il:3),surface(i)%sp04,ires)
   if (iproc==0) call scanfile (params%infile,'surface_param_05_'//cm(il:3),surface(i)%sp05,ires)
   if (iproc==0) call scanfile (params%infile,'surface_param_06_'//cm(il:3),surface(i)%sp06,ires)
   if (iproc==0) call scanfile (params%infile,'surface_param_07_'//cm(il:3),surface(i)%sp07,ires)
   if (iproc==0) call scanfile (params%infile,'surface_param_08_'//cm(il:3),surface(i)%sp08,ires)
   if (iproc==0) call scanfile (params%infile,'surface_param_09_'//cm(il:3),surface(i)%sp09,ires)
   if (iproc==0) call scanfile (params%infile,'surface_param_10_'//cm(il:3),surface(i)%sp10,ires)
   
   call mpi_bcast(surface(i)%sp01,1,mpi_double_precision,0,mpi_comm_world,ierr)
   call mpi_bcast(surface(i)%sp02,1,mpi_double_precision,0,mpi_comm_world,ierr)
   call mpi_bcast(surface(i)%sp03,1,mpi_double_precision,0,mpi_comm_world,ierr)
   call mpi_bcast(surface(i)%sp04,1,mpi_double_precision,0,mpi_comm_world,ierr)
   call mpi_bcast(surface(i)%sp05,1,mpi_double_precision,0,mpi_comm_world,ierr)
   call mpi_bcast(surface(i)%sp06,1,mpi_double_precision,0,mpi_comm_world,ierr)
   call mpi_bcast(surface(i)%sp07,1,mpi_double_precision,0,mpi_comm_world,ierr)
   call mpi_bcast(surface(i)%sp08,1,mpi_double_precision,0,mpi_comm_world,ierr)
   call mpi_bcast(surface(i)%sp09,1,mpi_double_precision,0,mpi_comm_world,ierr)
   call mpi_bcast(surface(i)%sp10,1,mpi_double_precision,0,mpi_comm_world,ierr)

   write(threadinfo%Logunit,*) 'surface ',i,'sp01',surface(i)%sp01
   write(threadinfo%Logunit,*) 'surface ',i,'sp02',surface(i)%sp02
   write(threadinfo%Logunit,*) 'surface ',i,'sp03',surface(i)%sp03
   write(threadinfo%Logunit,*) 'surface ',i,'sp04',surface(i)%sp04
   write(threadinfo%Logunit,*) 'surface ',i,'sp05',surface(i)%sp05
   write(threadinfo%Logunit,*) 'surface ',i,'sp06',surface(i)%sp06
   write(threadinfo%Logunit,*) 'surface ',i,'sp07',surface(i)%sp07
   write(threadinfo%Logunit,*) 'surface ',i,'sp08',surface(i)%sp08
   write(threadinfo%Logunit,*) 'surface ',i,'sp09',surface(i)%sp09
   write(threadinfo%Logunit,*) 'surface ',i,'sp10',surface(i)%sp10

   surface(i)%activation_time=-1.d0
   if (iproc==0) call scanfile (params%infile,'activation_time_'//cm(il:3),surface(i)%activation_time,ires)
   call mpi_bcast(surface(i)%activation_time,1,mpi_double_precision,0,mpi_comm_world,ierr)
   write(threadinfo%Logunit,*) 'surface ',i,'activation_time',surface(i)%activation_time

   surface(i)%leveloct=6
   if (iproc==0) call scanfile (params%infile,'leveloct'//cm(il:3),surface(i)%leveloct,ires)
   call mpi_bcast(surface(i)%leveloct,1,mpi_integer,0,mpi_comm_world,ierr)
   write(threadinfo%Logunit,*) 'surface ',i,'leveloct',surface(i)%leveloct

enddo

!=====[face refinement parameters]=============================================

params%ref_on_faces=.false.
if(iproc==0) then
   call scanfile (params%infile,'ref_on_faces',answer,ires)
   params%ref_on_faces=(trim(answer)=='T')
end if
call mpi_bcast(params%ref_on_faces,1,mpi_logical,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'ref_on_faces',params%ref_on_faces

do i=1,6
   write(cm,'(i3)') i

   if (iproc==0) call scanfile (params%infile,'level_face'//cm(3:3),cube_faces(i)%level,ires)
   call mpi_bcast(cube_faces(i)%level,1,mpi_integer,0,mpi_comm_world,ierr)
   write(threadinfo%Logunit,*) 'cube_faces level',cube_faces(i)%level

   if (iproc==0) call scanfile (params%infile,'l'//cm(3:3),cube_faces(i)%l,ires)
   call mpi_bcast(cube_faces(i)%l,1,mpi_double_precision,0,mpi_comm_world,ierr)
   write(threadinfo%Logunit,*) 'cube_faces l',cube_faces(i)%l

   if (iproc==0) call scanfile (params%infile,'r'//cm(3:3),cube_faces(i)%r,ires)
   call mpi_bcast(cube_faces(i)%r,1,mpi_double_precision,0,mpi_comm_world,ierr)
   write(threadinfo%Logunit,*) 'cube_faces r',cube_faces(i)%r

   if (iproc==0) call scanfile (params%infile,'b'//cm(3:3),cube_faces(i)%b,ires)
   call mpi_bcast(cube_faces(i)%b,1,mpi_double_precision,0,mpi_comm_world,ierr)
   write(threadinfo%Logunit,*) 'cube_faces b',cube_faces(i)%b

   if (iproc==0) call scanfile (params%infile,'t'//cm(3:3),cube_faces(i)%t,ires)
   call mpi_bcast(cube_faces(i)%t,1,mpi_double_precision,0,mpi_comm_world,ierr)
   write(threadinfo%Logunit,*) 'cube_faces t',cube_faces(i)%t
end do

params%noctreemax=100000
if(iproc==0) call scanfile (params%infile,'noctreemax',params%noctreemax,ires)
call mpi_bcast(params%noctreemax,1,mpi_integer,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'noctreemax',params%noctreemax

params%nonlinear_iterations=3
if(iproc==0) call scanfile (params%infile,'nonlinear_iterations',params%nonlinear_iterations,ires)
call mpi_bcast(params%nonlinear_iterations,1,mpi_integer,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'nonlinear_iterations',params%nonlinear_iterations

params%initial_refine_level=6
if (iproc==0) call scanfile (params%infile,'initial_refine_level',params%initial_refine_level,ires)
call mpi_bcast(params%initial_refine_level,1,mpi_integer,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'initial_refine_level',params%initial_refine_level

params%dt=0.5d0
if (iproc==0) call scanfile (params%infile,'dt',params%dt,ires)
call mpi_bcast(params%dt,1,mpi_double_precision,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'dt',params%dt

params%penalty=1.d8
if (iproc==0) call scanfile (params%infile,'penalty',params%penalty,ires)
call mpi_bcast(params%penalty,1,mpi_double_precision,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'penalty',params%penalty

params%tempscale=1.d0
if (iproc==0) call scanfile (params%infile,'tempscale',params%tempscale,ires)
call mpi_bcast(params%tempscale,1,mpi_double_precision,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'tempscale',params%tempscale
      
params%refine_ratio=1.d0                                                          
if (iproc==0) call scanfile (params%infile,'refine_ratio',params%refine_ratio,ires)              
call mpi_bcast(params%refine_ratio,1,mpi_double_precision,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'refine_ratio',params%refine_ratio

params%refine_criterion=1                                                       
if(iproc==0) call scanfile (params%infile,'refine_criterion',params%refine_criterion,ires)    
call mpi_bcast(params%refine_criterion,1,mpi_integer,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'refine_criterion',params%refine_criterion

params%octree_refine_ratio=1.d0                                                   
if (iproc.eq.0) call scanfile (params%infile,'octree_refine_ratio',params%octree_refine_ratio,ires)
call mpi_bcast(params%octree_refine_ratio,1,mpi_double_precision,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'octree_refine_ratio',params%octree_refine_ratio


params%courant=.5d0
if( iproc==0) call scanfile (params%infile,'courant',params%courant,ires)
call mpi_bcast(params%courant,1,mpi_double_precision,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'courant',params%courant

!//////////////////////////////////////

params%npmin=8
if (iproc==0) call scanfile (params%infile,'npmin',params%npmin,ires)
call mpi_bcast(params%npmin,1,mpi_integer,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'npmin',params%npmin

params%npmax=16
if (iproc==0) call scanfile (params%infile,'npmax',params%npmax,ires)
call mpi_bcast(params%npmax,1,mpi_integer,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'npmax',params%npmax

!//////////////////////////////////////

params%griditer=-10
if (iproc==0) call scanfile (params%infile,'griditer',params%griditer,ires)
call mpi_bcast(params%griditer,1,mpi_integer,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'griditer',params%griditer

params%tol=1.d-3
if (iproc==0) call scanfile (params%infile,'tol',params%tol,ires)
call mpi_bcast(params%tol,1,mpi_double_precision,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'tol',params%tol

!//////////////////////////////////////

!params%distance_exponent=-log(2.d0)/log(cos(params%anglemax))
params%distance_exponent=1.d0
write(threadinfo%Logunit,*) 'distance_exponent',params%distance_exponent

!//////////////////////////////////////

params%niter_move=10
if (iproc==0) call scanfile (params%infile,'niter_move',params%niter_move,ires)
call mpi_bcast(params%niter_move,1,mpi_integer,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'niter_move',params%niter_move

!//////////////////////////////////////

params%irestart=0
if (iproc==0) call scanfile (params%infile,'irestart',params%irestart,ires)
call mpi_bcast(params%irestart,1,mpi_integer,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'irestart',params%irestart

!//////////////////////////////////////

params%restartfile='OUT/xxxxxx.bin'
if(iproc==0) call scanfile (params%infile,'restartfile',params%restartfile,ires)
call mpi_bcast(params%restartfile,128,mpi_character,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'restartfile',params%restartfile

!//////////////////////////////////////

params%ismooth=.false.
if (iproc==0) then
   call scanfile (params%infile,'ismooth',answer,ires)
   params%ismooth = (trim(answer)=='T')
end if
call mpi_bcast(params%ismooth,1,mpi_logical,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'ismooth',params%ismooth

!//////////////////////////////////////

params%nboxes=0
if (iproc==0) call scanfile (params%infile,'nboxes',params%nboxes,ires)
call mpi_bcast(params%nboxes,1,mpi_integer,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'nboxes',params%nboxes

!//////////////////////////////////////

params%nb_iter_nl_min=0
if (iproc==0) call scanfile (params%infile,'nb_iter_nl_min',params%nb_iter_nl_min,ires)
call mpi_bcast(params%nb_iter_nl_min,1,mpi_integer,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'nb_iter_nl_min',params%nb_iter_nl_min

!//////////////////////////////////////

params%visualise_matrix=.false.
if (iproc==0) then
   call scanfile (params%infile,'visualise_matrix',answer,ires)
   params%visualise_matrix=(trim(answer)=='T')
end if
call mpi_bcast(params%visualise_matrix,1,mpi_logical,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'visualise_matrix',params%visualise_matrix

!//////////////////////////////////////

params%renumber_nodes  = .false.
if (iproc==0) then
   call scanfile (params%infile,'renumber_nodes',answer,ires)
   params%renumber_nodes=(trim(answer)=='T')
end if
call mpi_bcast(params%renumber_nodes,1,mpi_logical,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'renumber_nodes',params%renumber_nodes

!//////////////////////////////////////

do i=1,params%nboxes

   write(cm,'(i3)') i
   il=1
   if (i.lt.100) il=2
   if (il.lt.10) il=3

   boxes(i)%x0=0.d0
   if (iproc==0) call scanfile (params%infile,'box'//cm(il:3)//'x0',boxes(i)%x0,ires)
   call mpi_bcast(boxes(i)%x0,1,mpi_double_precision,0,mpi_comm_world,ierr)
   write(threadinfo%Logunit,*) 'box x0 ',boxes(i)%x0

   boxes(i)%x1=0.d0
   if (iproc==0) call scanfile (params%infile,'box'//cm(il:3)//'x1',boxes(i)%x1,ires)
   call mpi_bcast(boxes(i)%x1,1,mpi_double_precision,0,mpi_comm_world,ierr)
   write(threadinfo%Logunit,*) 'box x1 ',boxes(i)%x1

   boxes(i)%y0=0.d0
   if (iproc==0) call scanfile (params%infile,'box'//cm(il:3)//'y0',boxes(i)%y0,ires)
   call mpi_bcast(boxes(i)%y0,1,mpi_double_precision,0,mpi_comm_world,ierr)
   write(threadinfo%Logunit,*) 'box y0 ',boxes(i)%y0

   boxes(i)%y1=0.d0
   if (iproc==0) call scanfile (params%infile,'box'//cm(il:3)//'y1',boxes(i)%y1,ires)
   call mpi_bcast(boxes(i)%y1,1,mpi_double_precision,0,mpi_comm_world,ierr)
   write(threadinfo%Logunit,*) 'box y1 ',boxes(i)%y1

   boxes(i)%z0=0.d0
   if (iproc==0) call scanfile (params%infile,'box'//cm(il:3)//'z0',boxes(i)%z0,ires)
   call mpi_bcast(boxes(i)%z0,1,mpi_double_precision,0,mpi_comm_world,ierr)
   write(threadinfo%Logunit,*) 'box z0 ',boxes(i)%z0

   boxes(i)%z1=0.d0
   if (iproc==0) call scanfile (params%infile,'box'//cm(il:3)//'z1',boxes(i)%z1,ires)
   call mpi_bcast(boxes(i)%z1,1,mpi_double_precision,0,mpi_comm_world,ierr)
   write(threadinfo%Logunit,*) 'box z1 ',boxes(i)%z1

   boxes(i)%level=1
   if (iproc==0) call scanfile (params%infile,'box'//cm(il:3)//'level',boxes(i)%level,ires)
   call mpi_bcast(boxes(i)%level,1,mpi_integer,0,mpi_comm_world,ierr)
   write(threadinfo%Logunit,*) 'box level',boxes(i)%level

enddo

params%nsections=0
if (iproc==0) call scanfile (params%infile,'nsections',params%nsections,ires)
call mpi_bcast(params%nsections,1,mpi_integer,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'nsections',params%nsections

do i=1,params%nsections

   write(cm,'(i3)') i
   il=1
   if (i.lt.100) il=2
   if (il.lt.10) il=3

   sections(i)%xyz=0
   if (iproc==0) call scanfile (params%infile,'xyz_'//cm(il:3),sections(i)%xyz,ires)

   if (iproc==0) call scanfile (params%infile,'slice_'//cm(il:3),sections(i)%slice,ires)
!   if (sections(i)%slice<0.d0 .or. sections(i)%slice>1.d0) call stop_run ('pb with slice value')

   sections(i)%flag_press = .false. 
   if (iproc==0) then
      call scanfile (params%infile,'flag_press_'//cm(il:3),answer,ires)
      sections(i)%flag_press=(trim(answer)=='T')
   end if

   sections(i)%flag_spress = .false. 
   if (iproc==0) then
      call scanfile (params%infile,'flag_spress_'//cm(il:3),answer,ires)
      sections(i)%flag_spress=(trim(answer)=='T')
   end if

   sections(i)%flag_e2d   = .false.
   if (iproc==0) then
      call scanfile (params%infile,'flag_e2d_'//cm(il:3),answer,ires)
      sections(i)%flag_e2d=(trim(answer)=='T')
   end if

   sections(i)%flag_e3d   = .false.
   if (iproc==0) then
      call scanfile (params%infile,'flag_e3d_'//cm(il:3),answer,ires)
      sections(i)%flag_e3d=(trim(answer)=='T')
   end if

   sections(i)%flag_strain   = .false.
   if (iproc==0) then
      call scanfile (params%infile,'flag_strain_'//cm(il:3),answer,ires)
      sections(i)%flag_strain=(trim(answer)=='T')
   end if

   sections(i)%flag_crit  = .false.
   if (iproc==0) then
      call scanfile (params%infile,'flag_crit_'//cm(il:3),answer,ires)
      sections(i)%flag_crit=(trim(answer)=='T')
   end if

   sections(i)%flag_grid  = .false.
   if (iproc==0) then
      call scanfile (params%infile,'flag_grid_'//cm(il:3),answer,ires)
      sections(i)%flag_grid=(trim(answer)=='T')
   end if

   sections(i)%flag_u  = .false.
   if (iproc==0) then
      call scanfile (params%infile,'flag_u_'//cm(il:3),answer,ires)
      sections(i)%flag_u=(trim(answer)=='T')
   end if
 
   sections(i)%flag_v  = .false.
   if (iproc==0) then
      call scanfile (params%infile,'flag_v_'//cm(il:3),answer,ires)
      sections(i)%flag_v=(trim(answer)=='T')
   end if

   sections(i)%flag_w  = .false.
   if (iproc==0) then
      call scanfile (params%infile,'flag_w_'//cm(il:3),answer,ires)
      sections(i)%flag_w=(trim(answer)=='T')
   end if

   sections(i)%flag_uvw  = .false.
   if (iproc==0) then
      call scanfile (params%infile,'flag_uvw_'//cm(il:3),answer,ires)
      sections(i)%flag_uvw=(trim(answer)=='T')
   end if

   sections(i)%flag_colour  = .false.
   if (iproc==0) then
      call scanfile (params%infile,'flag_colour_'//cm(il:3),answer,ires)
      sections(i)%flag_colour=(trim(answer)=='T')
   end if

   sections(i)%flag_mu  = .false.
   if (iproc==0) then
      call scanfile (params%infile,'flag_mu_'//cm(il:3),answer,ires)
      sections(i)%flag_mu=(trim(answer)=='T')
   end if

   sections(i)%flag_plastic  = .false.
   if (iproc==0) then
      call scanfile (params%infile,'flag_plastic_'//cm(il:3),answer,ires)
      sections(i)%flag_plastic=(trim(answer)=='T')
   end if

   sections(i)%flag_q  = .false.
   if (iproc==0) then
      call scanfile (params%infile,'flag_q_'//cm(il:3),answer,ires)
      sections(i)%flag_q=(trim(answer)=='T')
   end if

   sections(i)%flag_lode  = .false.
   if (iproc==0) then
      call scanfile (params%infile,'flag_lode_'//cm(il:3),answer,ires)
      sections(i)%flag_lode=(trim(answer)=='T')
   end if

   sections(i)%flag_vfield  = .false.
   if (iproc==0) then
      call scanfile (params%infile,'flag_vfield_'//cm(il:3),answer,ires)
      sections(i)%flag_vfield=(trim(answer)=='T')
   end if

   sections(i)%flag_lsf  = .false.
   if (iproc==0) then
      call scanfile (params%infile,'flag_lsf_'//cm(il:3),answer,ires)
      sections(i)%flag_lsf=(trim(answer)=='T')
   end if

   sections(i)%flag_temp  = .false.
   if (iproc==0) then
      call scanfile (params%infile,'flag_temp_'//cm(il:3),answer,ires)
      sections(i)%flag_temp=(trim(answer)=='T')
   end if

   sections(i)%flag_velvect  = .false.
   if (iproc==0) then
      call scanfile (params%infile,'flag_velvect_'//cm(il:3),answer,ires)
      sections(i)%flag_velvect=(trim(answer)=='T')
   end if

   sections(i)%scale=100.d0
   if (iproc==0) call scanfile (params%infile,'scale_'//cm(il:3),sections(i)%scale,ires)

   if (iproc==0) call scanfile (params%infile,'colormap_'//cm(il:3),sections(i)%colormap,ires)

   if (iproc==0) call scanfile (params%infile,'ncolours_'//cm(il:3),sections(i)%ncolours,ires)

end do


!=====[erosion parameters]=====================================================

params%erosion=.false.
if (iproc==0) then
   call scanfile (params%infile,'erosion',answer,ires)
   params%erosion=(trim(answer)=='T')
end if
call mpi_bcast(params%erosion,1,mpi_logical,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'erosion',params%erosion

params%zerosion=0.d0
if (iproc==0) call scanfile (params%infile,'zerosion',params%zerosion,ires)
call mpi_bcast(params%zerosion,1,mpi_double_precision,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'zerosion',params%zerosion

params%length_scale=1.d0
if (iproc==0) call scanfile (params%infile,'length_scale',params%length_scale,ires)
call mpi_bcast(params%length_scale,1,mpi_double_precision,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'length_scale',params%length_scale

params%velocity_scale=1.d0
if (iproc==0) call scanfile (params%infile,'velocity_scale',params%velocity_scale,ires)
call mpi_bcast(params%velocity_scale,1,mpi_double_precision,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'velocity_scale',params%velocity_scale

params%fluvial_erosion=2.d-1
if (iproc==0) call scanfile (params%infile,'fluvial_erosion',params%fluvial_erosion,ires)
call mpi_bcast(params%fluvial_erosion,1,mpi_double_precision,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'fluvial_erosion',params%fluvial_erosion

params%diffusion_erosion=8.d0
if (iproc==0) call scanfile (params%infile,'diffusion_erosion',params%diffusion_erosion,ires)
call mpi_bcast(params%diffusion_erosion,1,mpi_double_precision,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'diffusion_erosion',params%diffusion_erosion

params%baselevelx0=1
if (iproc==0) call scanfile (params%infile,'baselevelx0',params%baselevelx0,ires)
call mpi_bcast(params%baselevelx0,1,mpi_integer,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'baselevelx0',params%baselevelx0

params%baselevelx1=1
if (iproc==0) call scanfile (params%infile,'baselevelx1',params%baselevelx1,ires)
call mpi_bcast(params%baselevelx1,1,mpi_integer,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'baselevelx1',params%baselevelx1

params%baselevely0=1
if (iproc==0) call scanfile (params%infile,'baselevely0',params%baselevely0,ires)
call mpi_bcast(params%baselevely0,1,mpi_integer,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'baselevely0',params%baselevely0

params%baselevely1=1
if (iproc==0) call scanfile (params%infile,'baselevely1',params%baselevely1,ires)
call mpi_bcast(params%baselevely1,1,mpi_integer,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'baselevely1',params%baselevely1

params%compute_qpgram=.false.
if (iproc.eq.0) then
   call scanfile (params%infile,'compute_qpgram',answer,ires)
   params%compute_qpgram=(trim(answer)=='T')
end if
call mpi_bcast(params%compute_qpgram,1,mpi_logical,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'compute_qpgram',params%compute_qpgram

!=====[flexure parameters]=====================================================

params%isostasy=.false.
if (iproc.eq.0) then
   call scanfile (params%infile,'isostasy',answer,ires)
   params%isostasy=(trim(answer)=='T')
endif

params%flexure=.false.
if (iproc.eq.0) then
   call scanfile (params%infile,'flexure',answer,ires)
   params%flexure=(trim(answer)=='T')
endif

params%elastic_plate_thickness=20.d3
if (iproc==0) call scanfile (params%infile,'elastic_plate_thickness',params%elastic_plate_thickness,ires)
call mpi_bcast(params%elastic_plate_thickness,1,mpi_double_precision,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'elastic_plate_thickness',params%elastic_plate_thickness

params%density_difference=3.d3
if (iproc==0) call scanfile (params%infile,'density_difference',params%density_difference,ires)
call mpi_bcast(params%density_difference,1,mpi_double_precision,0,mpi_comm_world,ierr)
write(threadinfo%Logunit,*) 'density_difference',params%density_difference

if(iproc.eq.0) call flush(8)
call flush(threadinfo%logunit)

!if (iproc.eq.0) print*,(i,mat(i)%plasticity_type,mat(i)%plasticity_parameters(:),i=1,params%nmat)

return
end
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|