read_input_file.f90

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

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

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

! subroutine where default values for all input parameters other than the
! controlling 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 - read in read_controlling_parameters
! 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 - read in read_controlling_parameters
! 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
! matrule is the flag that determines the rule used for assigning material
!    properties within the model volume
! 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 initial 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 - read in
!    read_controlling_parameters
! 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 constants
use definitions
!use mpi
use threads

implicit none

include 'mpif.h'

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
type (nest_info) :: nest
type (bc_definition) :: bcdef

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

integer ires,i,il,ierr,iproc,nproc,j
character(len=72) :: shift
character(len=3)  :: cm
character(len=1)  :: answer

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

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

shift=' '

!==============================================================================
!==============================================================================

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)

params%nstep_spinup=0
if (iproc==0) call scanfile (params%infile,'nstep_spinup',params%nstep_spinup,ires)
call mpi_bcast(params%nstep_spinup,1,mpi_integer,0,mpi_comm_world,ierr)

params%savstep=1
if (iproc==0) call scanfile (params%infile,'savstep',params%savstep,ires)
call mpi_bcast(params%savstep,1,mpi_integer,0,mpi_comm_world,ierr)

params%savoffset=0
if (iproc==0) call scanfile (params%infile,'savoffset',params%savoffset,ires)
call mpi_bcast(params%savoffset,1,mpi_integer,0,mpi_comm_world,ierr)

!=====[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)

params%bulkvisc=.false.
if (iproc==0) then
  call scanfile (params%infile,'bulkvisc',answer,ires)
  if (ires==1) params%bulkvisc=(trim(answer)=='T')
endif
call mpi_bcast(params%bulkvisc,1,mpi_logical,0,mpi_comm_world,ierr)

params%init_e2d=.false.
if (iproc==0) then
  call scanfile (params%infile,'init_e2d',answer,ires)
  if (ires==1) params%init_e2d=(trim(answer)=='T')
endif
call mpi_bcast(params%init_e2d,1,mpi_logical,0,mpi_comm_world,ierr)

params%e2d0=1.d0
if (iproc==0) call scanfile (params%infile,'e2d0',params%e2d0,ires)
call mpi_bcast(params%e2d0,1,mpi_double_precision,0,mpi_comm_world,ierr)

params%plastic_stress_min=-1.d0
if (iproc==0) call scanfile (params%infile,'plastic_stress_min',params%plastic_stress_min,ires)
call mpi_bcast(params%plastic_stress_min,1,mpi_double_precision,0,mpi_comm_world,ierr)

params%pressure0=0.d0
if (iproc==0) call scanfile (params%infile,'pressure0',params%pressure0,ires)
call mpi_bcast(params%pressure0,1,mpi_double_precision,0,mpi_comm_world,ierr)

params%dommat=-1
if (iproc==0) call scanfile (params%infile,'dommat',params%dommat,ires)
call mpi_bcast(params%dommat,1,mpi_integer,0,mpi_comm_world,ierr)

params%domvol=1.d-10
if (iproc==0) call scanfile (params%infile,'domvol',params%domvol,ires)
call mpi_bcast(params%domvol,1,mpi_double_precision,0,mpi_comm_world,ierr)

params%materials_on_cloud=.false.
if (iproc==0) then
  call scanfile (params%infile,'materials_on_cloud',answer,ires)
  if (ires==1) params%materials_on_cloud=(trim(answer)=='T')
endif
call mpi_bcast(params%materials_on_cloud,1,mpi_logical,0,mpi_comm_world,ierr)

if (params%materials_on_cloud) then
  params%move_cloud_at_midpoint=.true.
else
  params%move_cloud_at_midpoint=.false.
endif
if (iproc==0) then
  call scanfile (params%infile,'move_cloud_at_midpoint',answer,ires)
  if (ires==1) params%move_cloud_at_midpoint=(trim(answer)=='T')
endif
call mpi_bcast(params%move_cloud_at_midpoint,1,mpi_logical,0,mpi_comm_world,ierr)

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

   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)

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

   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)

   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)

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

   mat(i)%fvisc_weak_type='tot_strain'
   if (iproc==0) call scanfile (params%infile,'fvisc_weak_type'//cm(il:3),mat(i)%fvisc_weak_type,ires)
   call mpi_bcast(mat(i)%fvisc_weak_type,16,mpi_character,0,mpi_comm_world,ierr)

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

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

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

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

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

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

   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)

   mat(i)%plasticity_ss_type_coh='tot_strain'
   if (iproc==0) call scanfile (params%infile,'plasticity_ss_type_coh'//cm(il:3),mat(i)%plasticity_ss_type_coh,ires)
   call mpi_bcast(mat(i)%plasticity_ss_type_coh,16,mpi_character,0,mpi_comm_world,ierr)

   mat(i)%plasticity_ss_type_phi='tot_strain'
   if (iproc==0) call scanfile (params%infile,'plasticity_ss_type_phi'//cm(il:3),mat(i)%plasticity_ss_type_phi,ires)
   call mpi_bcast(mat(i)%plasticity_ss_type_phi,16,mpi_character,0,mpi_comm_world,ierr)
   
   mat(i)%compactible=.false.
   if (iproc==0) then
     call scanfile (params%infile,'compactible'//cm(il:3),answer,ires)
     if (ires==1) mat(i)%compactible=(trim(answer)=='T')
   endif
   call mpi_bcast(mat(i)%compactible,1,mpi_logical,0,mpi_comm_world,ierr)
   
   mat(i)%grain_density=1.d0
   if (i.eq.0) mat(i)%grain_density=0.d0
   if (iproc==0) call scanfile (params%infile,'grain_density'//cm(il:3),mat(i)%grain_density,ires)
   call mpi_bcast(mat(i)%grain_density,1,mpi_double_precision,0,mpi_comm_world,ierr)
   
   mat(i)%pore_fluid_density=1.d0
   if (i.eq.0) mat(i)%pore_fluid_density=0.d0
   if (iproc==0) call scanfile (params%infile,'pore_fluid_density'//cm(il:3),mat(i)%pore_fluid_density,ires)
   call mpi_bcast(mat(i)%pore_fluid_density,1,mpi_double_precision,0,mpi_comm_world,ierr)
   
   mat(i)%surface_porosity=1.d0
   if (i.eq.0) mat(i)%surface_porosity=0.d0
   if (iproc==0) call scanfile (params%infile,'surface_porosity'//cm(il:3),mat(i)%surface_porosity,ires)
   call mpi_bcast(mat(i)%surface_porosity,1,mpi_double_precision,0,mpi_comm_world,ierr)
   
   mat(i)%compaction_coef=1.d0
   if (i.eq.0) mat(i)%compaction_coef=0.d0
   if (iproc==0) call scanfile (params%infile,'compaction_coef'//cm(il:3),mat(i)%compaction_coef,ires)
   call mpi_bcast(mat(i)%compaction_coef,1,mpi_double_precision,0,mpi_comm_world,ierr)

   mat(i)%plasticity_parameters=0.d0
   mat(i)%plasticity_parameters(3)=-1.d0
   mat(i)%plasticity_parameters(6)=-1.d0
   mat(i)%plasticity_parameters(9)=-1.d0
   mat(i)%plasticity_parameters(12)=-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)
     if (iproc==0) call scanfile (params%infile,'plasticity_10th_param'//cm(il:3),mat(i)%plasticity_parameters(10),ires)
     if (iproc==0) call scanfile (params%infile,'plasticity_11th_param'//cm(il:3),mat(i)%plasticity_parameters(11),ires)
     if (iproc==0) call scanfile (params%infile,'plasticity_12th_param'//cm(il:3),mat(i)%plasticity_parameters(12),ires)
     if (iproc==0) call scanfile (params%infile,'plasticity_13th_param'//cm(il:3),mat(i)%plasticity_parameters(13),ires)
     if (iproc==0) call scanfile (params%infile,'plasticity_14th_param'//cm(il:3),mat(i)%plasticity_parameters(14),ires)
     call mpi_bcast(mat(i)%plasticity_parameters,14,mpi_double_precision,0,mpi_comm_world,ierr)
   endif

   mat(i)%mattrans=-1.d0
   if (iproc==0) call scanfile (params%infile,'mattrans_xmin'//cm(il:3),mat(i)%mattrans(1),ires)
   if (iproc==0) call scanfile (params%infile,'mattrans_xmax'//cm(il:3),mat(i)%mattrans(2),ires)
   if (iproc==0) call scanfile (params%infile,'mattrans_ymin'//cm(il:3),mat(i)%mattrans(3),ires)
   if (iproc==0) call scanfile (params%infile,'mattrans_ymax'//cm(il:3),mat(i)%mattrans(4),ires)
   if (iproc==0) call scanfile (params%infile,'mattrans_zmin'//cm(il:3),mat(i)%mattrans(5),ires)
   if (iproc==0) call scanfile (params%infile,'mattrans_zmax'//cm(il:3),mat(i)%mattrans(6),ires)
   call mpi_bcast(mat(i)%mattrans,6,mpi_double_precision,0,mpi_comm_world,ierr)
   
   mat(i)%transnum=i
   if (iproc==0) then
     call scanfile (params%infile,'transnum_xmin'//cm(il:3),mat(i)%transnum(1),ires)
     call scanfile (params%infile,'transnum_xmax'//cm(il:3),mat(i)%transnum(2),ires)
     call scanfile (params%infile,'transnum_ymin'//cm(il:3),mat(i)%transnum(3),ires)
     call scanfile (params%infile,'transnum_ymax'//cm(il:3),mat(i)%transnum(4),ires)
     call scanfile (params%infile,'transnum_zmin'//cm(il:3),mat(i)%transnum(5),ires)
     call scanfile (params%infile,'transnum_zmax'//cm(il:3),mat(i)%transnum(6),ires)
     if (maxval(mat(i)%transnum).gt.params%nmat) call stop_run ('Error: transnum value is greater than nmat in read_input_file$')
   endif
   call mpi_bcast(mat(i)%transnum,6,mpi_integer,0,mpi_comm_world,ierr)
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)

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)

params%vex=1.d0
if (iproc==0) call scanfile (params%infile,'vex',params%vex,ires)
call mpi_bcast(params%vex,1,mpi_double_precision,0,mpi_comm_world,ierr)

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)

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)

params%matrule=0
if (iproc==0) call scanfile (params%infile,'matrule',params%matrule,ires)
call mpi_bcast(params%matrule,1,mpi_integer,0,mpi_comm_world,ierr)

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)

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)

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

params%ztemp=1.d0
if (iproc==0) call scanfile (params%infile,'ztemp',params%ztemp,ires)
params%ztemp=params%ztemp/params%vex
call mpi_bcast(params%ztemp,1,mpi_double_precision,0,mpi_comm_world,ierr)

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

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

params%sstemp_viscosity_spinup=-1.d0
if (iproc==0) call scanfile (params%infile,'sstemp_viscosity_spinup',params%sstemp_viscosity_spinup,ires)
call mpi_bcast(params%sstemp_viscosity_spinup,1,mpi_double_precision,0,mpi_comm_world,ierr)

params%sstemp_penalty_spinup=-1.d0
if (iproc==0) call scanfile (params%infile,'sstemp_penalty_spinup',params%sstemp_penalty_spinup,ires)
call mpi_bcast(params%sstemp_penalty_spinup,1,mpi_double_precision,0,mpi_comm_world,ierr)

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

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

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)

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

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

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

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

!=====[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)

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

   surface(i)%material_main=1
   if (iproc==0) call scanfile (params%infile,'material'//cm(il:3),surface(i)%material_main,ires)
   call mpi_bcast(surface(i)%material_main,1,mpi_integer,0,mpi_comm_world,ierr)

   surface(i)%material_spinup=surface(i)%material_main
   if (iproc==0) call scanfile (params%infile,'material_spinup'//cm(il:3),surface(i)%material_spinup,ires)
   call mpi_bcast(surface(i)%material_spinup,1,mpi_integer,0,mpi_comm_world,ierr)

   surface(i)%material=surface(i)%material_main
   call mpi_bcast(surface(i)%material,1,mpi_integer,0,mpi_comm_world,ierr)

   surface(i)%rand=.false.
   if (iproc==0) then
      call scanfile (params%infile,'rand'//cm(il:3),answer,ires)
      if (ires==1) surface(i)%rand=(trim(answer)=='T')
   end if
   call mpi_bcast(surface(i)%rand,1,mpi_logical,0,mpi_comm_world,ierr)

   surface(i)%surf_for_mat_props=.true.
   if (iproc==0) then
      call scanfile (params%infile,'surf_for_mat_props'//cm(il:3),answer,ires)
      if (ires==1) surface(i)%surf_for_mat_props=(trim(answer)=='T')
   end if
   call mpi_bcast(surface(i)%surf_for_mat_props,1,mpi_logical,0,mpi_comm_world,ierr)

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

   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)

   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)

   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)

   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)

   surface(i)%fixed_surf_spinup=.false.
   if (iproc==0) then
      call scanfile (params%infile,'fixed_surf_spinup'//cm(il:3),answer,ires)
      if (ires==1) surface(i)%fixed_surf_spinup=(trim(answer)=='T')
   end if
   call mpi_bcast(surface(i)%fixed_surf_spinup,1,mpi_logical,0,mpi_comm_world,ierr)

   surface(i)%fixed_surf=.false.
   if (iproc==0) then
      call scanfile (params%infile,'fixed_surf'//cm(il:3),answer,ires)
      if (ires==1) surface(i)%fixed_surf=(trim(answer)=='T')
   end if
   call mpi_bcast(surface(i)%fixed_surf,1,mpi_logical,0,mpi_comm_world,ierr)

   surface(i)%surface_type = 0
   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)
 
   surface(i)%sp01 = 0.d0
   if (iproc==0) call scanfile (params%infile,'surface_param_01_'//cm(il:3),surface(i)%sp01,ires)
   surface(i)%sp02 = 0.d0
   if (iproc==0) call scanfile (params%infile,'surface_param_02_'//cm(il:3),surface(i)%sp02,ires)
   surface(i)%sp03 = 0.d0
   if (iproc==0) call scanfile (params%infile,'surface_param_03_'//cm(il:3),surface(i)%sp03,ires)
   surface(i)%sp04 = 0.d0
   if (iproc==0) call scanfile (params%infile,'surface_param_04_'//cm(il:3),surface(i)%sp04,ires)
   surface(i)%sp05 = 0.d0
   if (iproc==0) call scanfile (params%infile,'surface_param_05_'//cm(il:3),surface(i)%sp05,ires)
   surface(i)%sp06 = 0.d0
   if (iproc==0) call scanfile (params%infile,'surface_param_06_'//cm(il:3),surface(i)%sp06,ires)
   surface(i)%sp07 = 0.d0
   if (iproc==0) call scanfile (params%infile,'surface_param_07_'//cm(il:3),surface(i)%sp07,ires)
   surface(i)%sp08 = 0.d0
   if (iproc==0) call scanfile (params%infile,'surface_param_08_'//cm(il:3),surface(i)%sp08,ires)
   surface(i)%sp09 = 0.d0
   if (iproc==0) call scanfile (params%infile,'surface_param_09_'//cm(il:3),surface(i)%sp09,ires)
   surface(i)%sp10 = 0.d0
   if (iproc==0) call scanfile (params%infile,'surface_param_10_'//cm(il:3),surface(i)%sp10,ires)
   surface(i)%sp11 = surface(i)%sp01
   if (iproc==0) call scanfile (params%infile,'surface_param_11_'//cm(il:3),surface(i)%sp11,ires)
   surface(i)%sp12 = 45.d0
   if (iproc==0) call scanfile (params%infile,'surface_param_12_'//cm(il:3),surface(i)%sp12,ires)
   surface(i)%sp13 = surface(i)%sp11
   if (iproc==0) call scanfile (params%infile,'surface_param_13_'//cm(il:3),surface(i)%sp13,ires)
   surface(i)%sp14 = 45.d0
   if (iproc==0) call scanfile (params%infile,'surface_param_14_'//cm(il:3),surface(i)%sp14,ires)
   surface(i)%sp15 = 1.d0
   if (iproc==0) call scanfile (params%infile,'surface_param_15_'//cm(il:3),surface(i)%sp15,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)
   call mpi_bcast(surface(i)%sp11,1,mpi_double_precision,0,mpi_comm_world,ierr)
   call mpi_bcast(surface(i)%sp12,1,mpi_double_precision,0,mpi_comm_world,ierr)
   call mpi_bcast(surface(i)%sp13,1,mpi_double_precision,0,mpi_comm_world,ierr)
   call mpi_bcast(surface(i)%sp14,1,mpi_double_precision,0,mpi_comm_world,ierr)
   call mpi_bcast(surface(i)%sp15,1,mpi_double_precision,0,mpi_comm_world,ierr)

   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)

   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)

   surface(i)%remove_after_mat_def=.false.
   if (iproc==0) then
      call scanfile (params%infile,'remove_after_mat_def'//cm(il:3),answer,ires)
      if (ires==1) surface(i)%remove_after_mat_def=(trim(answer)=='T')
   end if
   call mpi_bcast(surface(i)%remove_after_mat_def,1,mpi_logical,0,mpi_comm_world,ierr)
enddo

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

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

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

   cube_faces(i)%level = 1
   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)

   cube_faces(i)%l = 0.d0
   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)

   cube_faces(i)%r = 0.d0
   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)

   cube_faces(i)%b = 0.d0
   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)

   cube_faces(i)%t = 0.d0
   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)
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)

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

params%nonlinear_iterations_spinup=params%nonlinear_iterations_main
if(iproc==0) call scanfile (params%infile,'nonlinear_iterations_spinup',params%nonlinear_iterations_spinup,ires)
call mpi_bcast(params%nonlinear_iterations_spinup,1,mpi_integer,0,mpi_comm_world,ierr)

params%nonlinear_iterations=params%nonlinear_iterations_main
call mpi_bcast(params%nonlinear_iterations,1,mpi_integer,0,mpi_comm_world,ierr)

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)

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

params%dt_spinup=params%dt_main
if (iproc==0) call scanfile (params%infile,'dt_spinup',params%dt_spinup,ires)
call mpi_bcast(params%dt_spinup,1,mpi_double_precision,0,mpi_comm_world,ierr)

params%dt=params%dt_main
call mpi_bcast(params%dt,1,mpi_double_precision,0,mpi_comm_world,ierr)

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

params%damp_factor=1.0d0
if (iproc==0) call scanfile (params%infile,'damp_factor',params%damp_factor,ires)
call mpi_bcast(params%damp_factor,1,mpi_double_precision,0,mpi_comm_world,ierr)

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

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

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)

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)

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)

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)

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)

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)

!=====[boundary conditions]=====================================================

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

bcdef%bctype=''
if (iproc.eq.0) call scanfile (params%infile,'bctype',bcdef%bctype,ires)
call mpi_bcast(bcdef%bctype,40,mpi_character,0,mpi_comm_world,ierr)

select case(trim(bcdef%bctype))
case('basic')
  bcdef%bcorder='xyz'
  if (iproc.eq.0) call scanfile (params%infile,'bcorder',bcdef%bcorder,ires)
  call mpi_bcast(bcdef%bcorder,3,mpi_character,0,mpi_comm_world,ierr)

  bcdef%fixux0=.true.
  bcdef%fixux1=.true.
  bcdef%fixvx0=.true.
  bcdef%fixvx1=.true.
  bcdef%fixwx0=.true.
  bcdef%fixwx1=.true.
  bcdef%fixuy0=.true.
  bcdef%fixuy1=.true.
  bcdef%fixvy0=.true.
  bcdef%fixvy1=.true.
  bcdef%fixwy0=.true.
  bcdef%fixwy1=.true.
  bcdef%fixuz0=.true.
  bcdef%fixuz1=.true.
  bcdef%fixvz0=.true.
  bcdef%fixvz1=.true.
  bcdef%fixwz0=.true.
  bcdef%fixwz1=.true.
  if(iproc==0) then
    call scanfile (params%infile,'fixux0',answer,ires)
    if (ires==1) bcdef%fixux0=(trim(answer)=='T')
  endif
  if(iproc==0) then
    call scanfile (params%infile,'fixux1',answer,ires)
    if (ires==1) bcdef%fixux1=(trim(answer)=='T')
  endif
  if(iproc==0) then
    call scanfile (params%infile,'fixvx0',answer,ires)
    if (ires==1) bcdef%fixvx0=(trim(answer)=='T')
  endif
  if(iproc==0) then
    call scanfile (params%infile,'fixvx1',answer,ires)
    if (ires==1) bcdef%fixvx1=(trim(answer)=='T')
  endif
  if(iproc==0) then
    call scanfile (params%infile,'fixwx0',answer,ires)
    if (ires==1) bcdef%fixwx0=(trim(answer)=='T')
  endif
  if(iproc==0) then
    call scanfile (params%infile,'fixwx1',answer,ires)
    if (ires==1) bcdef%fixwx1=(trim(answer)=='T')
  endif
  if(iproc==0) then
    call scanfile (params%infile,'fixuy0',answer,ires)
    if (ires==1) bcdef%fixuy0=(trim(answer)=='T')
  endif
  if(iproc==0) then
    call scanfile (params%infile,'fixuy1',answer,ires)
    if (ires==1) bcdef%fixuy1=(trim(answer)=='T')
  endif
  if(iproc==0) then
    call scanfile (params%infile,'fixvy0',answer,ires)
    if (ires==1) bcdef%fixvy0=(trim(answer)=='T')
  endif
  if(iproc==0) then
    call scanfile (params%infile,'fixvy1',answer,ires)
    if (ires==1) bcdef%fixvy1=(trim(answer)=='T')
  endif
  if(iproc==0) then
    call scanfile (params%infile,'fixwy0',answer,ires)
    if (ires==1) bcdef%fixwy0=(trim(answer)=='T')
  endif
  if(iproc==0) then
    call scanfile (params%infile,'fixwy1',answer,ires)
    if (ires==1) bcdef%fixwy1=(trim(answer)=='T')
  endif
  if(iproc==0) then
    call scanfile (params%infile,'fixuz0',answer,ires)
    if (ires==1) bcdef%fixuz0=(trim(answer)=='T')
  endif
  if(iproc==0) then
    call scanfile (params%infile,'fixuz1',answer,ires)
    if (ires==1) bcdef%fixuz1=(trim(answer)=='T')
  endif
  if(iproc==0) then
    call scanfile (params%infile,'fixvz0',answer,ires)
    if (ires==1) bcdef%fixvz0=(trim(answer)=='T')
  endif
  if(iproc==0) then
    call scanfile (params%infile,'fixvz1',answer,ires)
    if (ires==1) bcdef%fixvz1=(trim(answer)=='T')
  endif
  if(iproc==0) then
    call scanfile (params%infile,'fixwz0',answer,ires)
    if (ires==1) bcdef%fixwz0=(trim(answer)=='T')
  endif
  if(iproc==0) then
    call scanfile (params%infile,'fixwz1',answer,ires)
    if (ires==1) bcdef%fixwz1=(trim(answer)=='T')
  endif
  call mpi_bcast(bcdef%fixux0,1,mpi_logical,0,mpi_comm_world,ierr)
  call mpi_bcast(bcdef%fixux1,1,mpi_logical,0,mpi_comm_world,ierr)
  call mpi_bcast(bcdef%fixvx0,1,mpi_logical,0,mpi_comm_world,ierr)
  call mpi_bcast(bcdef%fixvx1,1,mpi_logical,0,mpi_comm_world,ierr)
  call mpi_bcast(bcdef%fixwx0,1,mpi_logical,0,mpi_comm_world,ierr)
  call mpi_bcast(bcdef%fixwx1,1,mpi_logical,0,mpi_comm_world,ierr)
  call mpi_bcast(bcdef%fixuy0,1,mpi_logical,0,mpi_comm_world,ierr)
  call mpi_bcast(bcdef%fixuy1,1,mpi_logical,0,mpi_comm_world,ierr)
  call mpi_bcast(bcdef%fixvy0,1,mpi_logical,0,mpi_comm_world,ierr)
  call mpi_bcast(bcdef%fixvy1,1,mpi_logical,0,mpi_comm_world,ierr)
  call mpi_bcast(bcdef%fixwy0,1,mpi_logical,0,mpi_comm_world,ierr)
  call mpi_bcast(bcdef%fixwy1,1,mpi_logical,0,mpi_comm_world,ierr)
  call mpi_bcast(bcdef%fixuz0,1,mpi_logical,0,mpi_comm_world,ierr)
  call mpi_bcast(bcdef%fixuz1,1,mpi_logical,0,mpi_comm_world,ierr)
  call mpi_bcast(bcdef%fixvz0,1,mpi_logical,0,mpi_comm_world,ierr)
  call mpi_bcast(bcdef%fixvz1,1,mpi_logical,0,mpi_comm_world,ierr)
  call mpi_bcast(bcdef%fixwz0,1,mpi_logical,0,mpi_comm_world,ierr)
  call mpi_bcast(bcdef%fixwz1,1,mpi_logical,0,mpi_comm_world,ierr)

  bcdef%ux0=0.d0
  bcdef%ux1=0.d0
  bcdef%vx0=0.d0
  bcdef%vx1=0.d0
  bcdef%wx0=0.d0
  bcdef%wx1=0.d0
  bcdef%uy0=0.d0
  bcdef%uy1=0.d0
  bcdef%vy0=0.d0
  bcdef%vy1=0.d0
  bcdef%wy0=0.d0
  bcdef%wy1=0.d0
  bcdef%uz0=0.d0
  bcdef%uz1=0.d0
  bcdef%vz0=0.d0
  bcdef%vz1=0.d0
  bcdef%wz0=0.d0
  bcdef%wz1=0.d0
  if (iproc==0) call scanfile (params%infile,'ux0',bcdef%ux0,ires)
  if (iproc==0) call scanfile (params%infile,'ux1',bcdef%ux1,ires)
  if (iproc==0) call scanfile (params%infile,'vx0',bcdef%vx0,ires)
  if (iproc==0) call scanfile (params%infile,'vx1',bcdef%vx1,ires)
  if (iproc==0) call scanfile (params%infile,'wx0',bcdef%wx0,ires)
  if (iproc==0) call scanfile (params%infile,'wx1',bcdef%wx1,ires)
  if (iproc==0) call scanfile (params%infile,'uy0',bcdef%uy0,ires)
  if (iproc==0) call scanfile (params%infile,'uy1',bcdef%uy1,ires)
  if (iproc==0) call scanfile (params%infile,'vy0',bcdef%vy0,ires)
  if (iproc==0) call scanfile (params%infile,'vy1',bcdef%vy1,ires)
  if (iproc==0) call scanfile (params%infile,'wy0',bcdef%wy0,ires)
  if (iproc==0) call scanfile (params%infile,'wy1',bcdef%wy1,ires)
  if (iproc==0) call scanfile (params%infile,'uz0',bcdef%uz0,ires)
  if (iproc==0) call scanfile (params%infile,'uz1',bcdef%uz1,ires)
  if (iproc==0) call scanfile (params%infile,'vz0',bcdef%vz0,ires)
  if (iproc==0) call scanfile (params%infile,'vz1',bcdef%vz1,ires)
  if (iproc==0) call scanfile (params%infile,'wz0',bcdef%wz0,ires)
  if (iproc==0) call scanfile (params%infile,'wz1',bcdef%wz1,ires)
  call mpi_bcast(bcdef%ux0,1,mpi_double_precision,0,mpi_comm_world,ierr)
  call mpi_bcast(bcdef%ux1,1,mpi_double_precision,0,mpi_comm_world,ierr)
  call mpi_bcast(bcdef%vx0,1,mpi_double_precision,0,mpi_comm_world,ierr)
  call mpi_bcast(bcdef%vx1,1,mpi_double_precision,0,mpi_comm_world,ierr)
  call mpi_bcast(bcdef%wx0,1,mpi_double_precision,0,mpi_comm_world,ierr)
  call mpi_bcast(bcdef%wx1,1,mpi_double_precision,0,mpi_comm_world,ierr)
  call mpi_bcast(bcdef%uy0,1,mpi_double_precision,0,mpi_comm_world,ierr)
  call mpi_bcast(bcdef%uy1,1,mpi_double_precision,0,mpi_comm_world,ierr)
  call mpi_bcast(bcdef%vy0,1,mpi_double_precision,0,mpi_comm_world,ierr)
  call mpi_bcast(bcdef%vy1,1,mpi_double_precision,0,mpi_comm_world,ierr)
  call mpi_bcast(bcdef%wy0,1,mpi_double_precision,0,mpi_comm_world,ierr)
  call mpi_bcast(bcdef%wy1,1,mpi_double_precision,0,mpi_comm_world,ierr)
  call mpi_bcast(bcdef%uz0,1,mpi_double_precision,0,mpi_comm_world,ierr)
  call mpi_bcast(bcdef%uz1,1,mpi_double_precision,0,mpi_comm_world,ierr)
  call mpi_bcast(bcdef%vz0,1,mpi_double_precision,0,mpi_comm_world,ierr)
  call mpi_bcast(bcdef%vz1,1,mpi_double_precision,0,mpi_comm_world,ierr)
  call mpi_bcast(bcdef%wz0,1,mpi_double_precision,0,mpi_comm_world,ierr)
  call mpi_bcast(bcdef%wz1,1,mpi_double_precision,0,mpi_comm_world,ierr)

case('segmented_s_line')
  bcdef%bc_parameters=0.d0
  if (iproc==0) call scanfile (params%infile,'bc_param1',bcdef%bc_parameters(1),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param2',bcdef%bc_parameters(2),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param3',bcdef%bc_parameters(3),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param4',bcdef%bc_parameters(4),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param5',bcdef%bc_parameters(5),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param6',bcdef%bc_parameters(6),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param7',bcdef%bc_parameters(7),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param8',bcdef%bc_parameters(8),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param9',bcdef%bc_parameters(9),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param10',bcdef%bc_parameters(10),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param11',bcdef%bc_parameters(11),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param12',bcdef%bc_parameters(12),ires)
  call mpi_bcast(bcdef%bc_parameters,14,mpi_double_precision,0,mpi_comm_world,ierr)

case('segmented_s_line_round')
  bcdef%bc_parameters=0.d0
  if (iproc==0) call scanfile (params%infile,'bc_param1',bcdef%bc_parameters(1),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param2',bcdef%bc_parameters(2),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param3',bcdef%bc_parameters(3),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param4',bcdef%bc_parameters(4),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param5',bcdef%bc_parameters(5),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param6',bcdef%bc_parameters(6),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param7',bcdef%bc_parameters(7),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param8',bcdef%bc_parameters(8),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param9',bcdef%bc_parameters(9),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param10',bcdef%bc_parameters(10),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param11',bcdef%bc_parameters(11),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param12',bcdef%bc_parameters(12),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param13',bcdef%bc_parameters(13),ires)
  call mpi_bcast(bcdef%bc_parameters,14,mpi_double_precision,0,mpi_comm_world,ierr)

case('rot_subduction')
  bcdef%bc_parameters=0.d0
  if (iproc==0) call scanfile (params%infile,'bc_param1',bcdef%bc_parameters(1),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param2',bcdef%bc_parameters(2),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param3',bcdef%bc_parameters(3),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param4',bcdef%bc_parameters(4),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param5',bcdef%bc_parameters(5),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param6',bcdef%bc_parameters(6),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param7',bcdef%bc_parameters(7),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param8',bcdef%bc_parameters(8),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param9',bcdef%bc_parameters(9),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param10',bcdef%bc_parameters(10),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param11',bcdef%bc_parameters(11),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param12',bcdef%bc_parameters(12),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param13',bcdef%bc_parameters(13),ires)
  call mpi_bcast(bcdef%bc_parameters,14,mpi_double_precision,0,mpi_comm_world,ierr)

case('rot_sub_init')
  bcdef%bc_parameters=0.d0
  if (iproc==0) call scanfile (params%infile,'bc_param1',bcdef%bc_parameters(1),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param2',bcdef%bc_parameters(2),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param3',bcdef%bc_parameters(3),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param4',bcdef%bc_parameters(4),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param5',bcdef%bc_parameters(5),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param6',bcdef%bc_parameters(6),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param7',bcdef%bc_parameters(7),ires)
  if (iproc==0) call scanfile (params%infile,'bc_param8',bcdef%bc_parameters(8),ires)