read_input_file.f90

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)

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)

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)

!=====[flexure parameters]=====================================================
params%isostasy=.false.
if (iproc.eq.0) then
   call scanfile (params%infile,'isostasy',answer,ires)
   params%isostasy=(trim(answer)=='T')
endif
call mpi_bcast(params%isostasy,1,mpi_logical,0,mpi_comm_world,ierr)

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

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

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)

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)

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

if (params%nest) then
   nest%lsoutfile='OUT/time_0001.bin'
   if (iproc.eq.0) call scanfile (params%infile,'lsoutfile',nest%lsoutfile,ires)
   call mpi_bcast(nest%lsoutfile,128,mpi_character,0,mpi_comm_world,ierr)

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

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

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

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

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

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

! Defined, but not broadcast or read from input file
!params%distance_exponent=-log(2.d0)/log(cos(params%anglemax))
params%distance_exponent=1.d0

! Input values are now written to stdout or log files here
if (params%debug.gt.0 .and. iproc.eq.0) then
  write(*,'(a)')         shift//'Input file values: '
  write(*,'(a,i4)')    shift//'nstep ',params%nstep
  write(*,'(a,i4)')    shift//'material0 ',material0
  write(*,'(a,l1)')    shift//'bulkvisc ',params%bulkvisc
  write(*,'(a,l1)')    shift//'init_e2d ',params%init_e2d
  write(*,'(a,e11.4)') shift//'e2d0 ',params%e2d0
  write(*,'(a,e11.4)') shift//'minp ',params%minp
  write(*,'(a,e11.4)') shift//'pressure0 ',params%pressure0
  do i=1,params%nmat
    write(*,'(a,i4,a)')  shift//'--- Properties for material ',i,' ---'
    write(*,'(a,e11.4)') shift//'density ',mat(i)%density
    write(*,'(a,e11.4)') shift//'viscosity ',mat(i)%viscosity
    write(*,'(a,e11.4)') shift//'penalty ',mat(i)%penalty
    write(*,'(a,e11.4)') shift//'expon ',mat(i)%expon
    write(*,'(a,e11.4)') shift//'activationenergy ',mat(i)%activationenergy
    write(*,'(a,e11.4)') shift//'expansion ',mat(i)%expansion
    write(*,'(a,e11.4)') shift//'diffusivity ',mat(i)%diffusivity
    write(*,'(a,e11.4)') shift//'heat ',mat(i)%heat
    write(*,'(a,a8)')    shift//'plasticity_type ',mat(i)%plasticity_type
    if (trim(mat(i)%plasticity_type).ne.'No') then
      write(*,'(a,i4)') shift//'Plasticity parameters for material ',i
      do j=1,9
        write(*,'(a,i1,a2,e11.4)') shift//'plasticity parameter ',j,' ',mat(i)%plasticity_parameters(j)
      enddo
    endif
    write (*,'(a,i4)') shift//'Material transitions for material ',i
    do j=1,6
      write(*,'(a,i1,a2,e11.4)') shift//'mattrans ',j,' ',mat(i)%mattrans(j)
    enddo
    do j=1,6
      write(*,'(a,i1,a2,i4)') shift//'transnum ',j,' ',mat(i)%transnum(j)
    enddo
  enddo
  write(*,'(a,e11.4)') shift//'viscositymin ',params%viscositymin
  write(*,'(a,e11.4)') shift//'viscositymax ',params%viscositymax
  write(*,'(a,e11.4)') shift//'vex ',params%vex
  write(*,'(a,i4)')    shift//'leveluniform_oct ',params%leveluniform_oct
  write(*,'(a,i4)')    shift//'levelmax_oct ',params%levelmax_oct
  write(*,'(a,i4)')    shift//'matrule ',params%matrule
  write(*,'(a,i4)')    shift//'levelcut ',params%levelcut
  write(*,'(a,i4)')    shift//'levelapprox ',params%levelapprox
  write(*,'(a,l1)')    shift//'calculate_temp ',params%calculate_temp
  write(*,'(a,e11.4)') shift//'ztemp ',params%ztemp
  write(*,'(a,i4)')    shift//'smoothing_type ',params%smoothing_type
  write(*,'(a,l1)')    shift//'normaladvect ',params%normaladvect
  write(*,'(a,l1)')    shift//'excl_vol ',params%excl_vol
  write(*,'(a,l1)')    shift//'adaptive_tol ',params%adaptive_tol
  do i=1,params%ns
    write(*,'(a,i4,a)')       shift//'--- Properties for surface ',i,' ---'
    write(*,'(a,i4,a,i4)')    shift//'surface ',i,' itype ',surface(i)%itype
    write(*,'(a,i4,a,i4)')    shift//'surface ',i,' material ',surface(i)%material
    write(*,'(a,i4,a,l1)')    shift//'surface ',i,' rand ',surface(i)%rand
    write(*,'(a,i4,a,i4)')    shift//'surface ',i,' levelt ',surface(i)%levelt
    write(*,'(a,i4,a,e11.4)') shift//'surface ',i,' stretch ',surface(i)%stretch
    write(*,'(a,i4,a,i4)')    shift//'surface ',i,' criterion ',surface(i)%criterion
    write(*,'(a,i4,a,e11.4)') shift//'surface ',i,' anglemax ',surface(i)%anglemax
    write(*,'(a,i4,a,e11.4)') shift//'surface ',i,' anglemaxoctree ',surface(i)%anglemaxoctree
    write(*,'(a,i4,a,i4)')    shift//'surface ',i,' spread_surface_points ',surface(i)%spread_surface_points
    write(*,'(a,i4,a,i4)')    shift//'surface ',i,' type ',surface(i)%surface_type
    write(*,'(a,i4,a,e11.4)') shift//'surface ',i,' sp01 ',surface(i)%sp01
    write(*,'(a,i4,a,e11.4)') shift//'surface ',i,' sp02 ',surface(i)%sp02
    write(*,'(a,i4,a,e11.4)') shift//'surface ',i,' sp03 ',surface(i)%sp03
    write(*,'(a,i4,a,e11.4)') shift//'surface ',i,' sp04 ',surface(i)%sp04
    write(*,'(a,i4,a,e11.4)') shift//'surface ',i,' sp05 ',surface(i)%sp05
    write(*,'(a,i4,a,e11.4)') shift//'surface ',i,' sp06 ',surface(i)%sp06
    write(*,'(a,i4,a,e11.4)') shift//'surface ',i,' sp07 ',surface(i)%sp07
    write(*,'(a,i4,a,e11.4)') shift//'surface ',i,' sp08 ',surface(i)%sp08
    write(*,'(a,i4,a,e11.4)') shift//'surface ',i,' sp09 ',surface(i)%sp09
    write(*,'(a,i4,a,e11.4)') shift//'surface ',i,' sp10 ',surface(i)%sp10
    write(*,'(a,i4,a,e11.4)') shift//'surface ',i,' sp11 ',surface(i)%sp11
    write(*,'(a,i4,a,e11.4)') shift//'surface ',i,' sp12 ',surface(i)%sp12
    write(*,'(a,i4,a,e11.4)') shift//'surface ',i,' activation_time ',surface(i)%activation_time
    write(*,'(a,i4,a,i4)')    shift//'surface ',i,' leveloct ',surface(i)%leveloct
  enddo
  write(*,'(a,l1)') shift//'ref_on_faces ',params%ref_on_faces
  if (params%ref_on_faces) then
    do i=1,6
      write(*,'(a,i1,a)')  shift//'--- Refinement values for face ',i,' ---'
      write(*,'(a,i4)')    shift//'cube_faces level ',cube_faces(i)%level
      write(*,'(a,e11.4)') shift//'cube_faces l ',cube_faces(i)%l
      write(*,'(a,e11.4)') shift//'cube_faces r ',cube_faces(i)%r
      write(*,'(a,e11.4)') shift//'cube_faces b ',cube_faces(i)%b
      write(*,'(a,e11.4)') shift//'cube_faces t ',cube_faces(i)%t
    enddo
  endif
  write(*,'(a,i12)') shift//'noctreemax ',params%noctreemax
  write(*,'(a,i4)') shift//'nonlinear_iterations ',params%nonlinear_iterations
  write(*,'(a,i4)') shift//'initial_refine_level ',params%initial_refine_level
  write(*,'(a,e11.4)') shift//'dt ',params%dt
  write(*,'(a,l1)') shift//'damp_surface ',params%damp_surface
  write(*,'(a,e11.4)') shift//'damp_factor ',params%damp_factor
  write(*,'(a,l1)') shift//'remove_surf_pts ',params%remove_surf_pts
  write(*,'(a,e11.4)') shift//'penalty ',params%penalty
  write(*,'(a,e11.4)') shift//'tempscale ',params%tempscale
  write(*,'(a,e11.4)') shift//'refine_ratio ',params%refine_ratio
  write(*,'(a,i4)') shift//'refine_criterion ',params%refine_criterion
  write(*,'(a,e11.4)') shift//'octree_refine_ratio ',params%octree_refine_ratio
  write(*,'(a,e11.4)') shift//'courant ',params%courant
  write(*,'(a,i4)') shift//'npmin ',params%npmin
  write(*,'(a,i4)') shift//'npmax ',params%npmax
  write(*,'(a,i4)') shift//'griditer ',params%griditer
  write(*,'(a,e11.4)') shift//'tol ',params%tol
  write(*,'(a,i4)') shift//'niter_move ',params%niter_move
  write(*,'(a,l1)') shift//'ismooth ',params%ismooth
  write(*,'(a,i4)') shift//'nb_iter_nl_min ',params%nb_iter_nl_min
  write(*,'(a,l1)') shift//'visualise_matrix ',params%visualise_matrix
  write(*,'(a,l1)') shift//'renumber_nodes ',params%renumber_nodes
  if (params%nboxes.gt.0) then
    do i=1,params%nboxes
      write(*,'(a,i3,a)')  shift//'--- Input values for box ',i,' ---'
      write(*,'(a,e11.4)') shift//'box x0 ',boxes(i)%x0
      write(*,'(a,e11.4)') shift//'box x1 ',boxes(i)%x1
      write(*,'(a,e11.4)') shift//'box y0 ',boxes(i)%y0
      write(*,'(a,e11.4)') shift//'box y1 ',boxes(i)%y1
      write(*,'(a,e11.4)') shift//'box z0 ',boxes(i)%z0
      write(*,'(a,e11.4)') shift//'box z1 ',boxes(i)%z1
      write(*,'(a,i4)')    shift//'box level ',boxes(i)%level
    enddo
  endif
  if (params%nsections.gt.0) then
    do i=1,params%nsections
      write(*,'(a,i3,a)')  shift//'--- Input values for section ',i,' ---'    
      write(*,'(a,i4)')    shift//'section xyz ',sections(i)%xyz
      write(*,'(a,e11.4)') shift//'section slice ',sections(i)%slice
      write(*,'(a,l1)')    shift//'section flag press ',sections(i)%flag_press
      write(*,'(a,l1)')    shift//'section flag spress ',sections(i)%flag_spress
      write(*,'(a,l1)')    shift//'section flag e2d ',sections(i)%flag_e2d
      write(*,'(a,l1)')    shift//'section flag e3d ',sections(i)%flag_e3d
      write(*,'(a,l1)')    shift//'section flag strain ',sections(i)%flag_strain
      write(*,'(a,l1)')    shift//'section flag crit ',sections(i)%flag_crit
      write(*,'(a,l1)')    shift//'section flag grid ',sections(i)%flag_grid
      write(*,'(a,l1)')    shift//'section flag u ',sections(i)%flag_u
      write(*,'(a,l1)')    shift//'section flag v ',sections(i)%flag_v  
      write(*,'(a,l1)')    shift//'section flag w ',sections(i)%flag_w 
      write(*,'(a,l1)')    shift//'section flag uvw ',sections(i)%flag_uvw
      write(*,'(a,l1)')    shift//'section flag colour ',sections(i)%flag_colour
      write(*,'(a,l1)')    shift//'section flag mu ',sections(i)%flag_mu
      write(*,'(a,l1)')    shift//'section flag plastic ',sections(i)%flag_plastic
      write(*,'(a,l1)')    shift//'section flag q ',sections(i)%flag_q
      write(*,'(a,l1)')    shift//'section flag lode ',sections(i)%flag_lode
      write(*,'(a,l1)')    shift//'section flag vfield ',sections(i)%flag_vfield
      write(*,'(a,l1)')    shift//'section flag lsf ',sections(i)%flag_lsf
      write(*,'(a,i3)')    shift//'section flag temp ',sections(i)%flag_temp
      write(*,'(a,l1)')    shift//'section flag velvect ',sections(i)%flag_velvect
      write(*,'(a,e11.4)') shift//'section scale ',sections(i)%scale
      write(*,'(a,a3)')    shift//'section colormap ',sections(i)%colormap
      write(*,'(a,i4)')    shift//'section ncolours ',sections(i)%ncolours
    enddo
  endif
  write(*,'(a,l1)') shift//'erosion ',params%erosion
  write(*,'(a,e11.4)') shift//'zerosion ',params%zerosion
  write(*,'(a,e11.4)') shift//'length_scale ',params%length_scale
  write(*,'(a,e11.4)') shift//'velocity_scale ',params%velocity_scale
  write(*,'(a,e11.4)') shift//'fluvial_erosion ',params%fluvial_erosion
  write(*,'(a,e11.4)') shift//'diffusion_erosion ',params%diffusion_erosion
  write(*,'(a,i4)') shift//'baselevelx0 ',params%baselevelx0
  write(*,'(a,i4)') shift//'baselevelx1 ',params%baselevelx1
  write(*,'(a,i4)') shift//'baselevely0 ',params%baselevely0
  write(*,'(a,i4)') shift//'baselevely1 ',params%baselevely1
  write(*,'(a,l1)') shift//'compute_qpgram ',params%compute_qpgram
  write(*,'(a)') shift//'--- Boundary condition parameters ---'
  write(*,'(a,l1)') shift//'invariants_2d ',params%invariants_2d
  write(*,'(a,a)') shift//'bctype ',bcdef%bctype
  select case(trim(bcdef%bctype))
  case('basic')
    write(*,'(a,a)') shift//'bcorder ',bcdef%bcorder
    write(*,'(a,l1)') shift//'fixux0 ',bcdef%fixux0
    write(*,'(a,l1)') shift//'fixux1 ',bcdef%fixux1
    write(*,'(a,l1)') shift//'fixvx0 ',bcdef%fixvx0
    write(*,'(a,l1)') shift//'fixvx1 ',bcdef%fixvx1
    write(*,'(a,l1)') shift//'fixwx0 ',bcdef%fixwx0
    write(*,'(a,l1)') shift//'fixwx1 ',bcdef%fixwx1
    write(*,'(a,l1)') shift//'fixuy0 ',bcdef%fixuy0
    write(*,'(a,l1)') shift//'fixuy1 ',bcdef%fixuy1
    write(*,'(a,l1)') shift//'fixvy0 ',bcdef%fixvy0
    write(*,'(a,l1)') shift//'fixvy1 ',bcdef%fixvy1
    write(*,'(a,l1)') shift//'fixwy0 ',bcdef%fixwy0
    write(*,'(a,l1)') shift//'fixwy1 ',bcdef%fixwy1
    write(*,'(a,l1)') shift//'fixuz0 ',bcdef%fixuz0
    write(*,'(a,l1)') shift//'fixuz1 ',bcdef%fixuz1
    write(*,'(a,l1)') shift//'fixvz0 ',bcdef%fixvz0
    write(*,'(a,l1)') shift//'fixvz1 ',bcdef%fixvz1
    write(*,'(a,l1)') shift//'fixwz0 ',bcdef%fixwz0
    write(*,'(a,l1)') shift//'fixwz1 ',bcdef%fixwz1
    write(*,'(a,e11.4)') shift//'ux0 ',bcdef%ux0
    write(*,'(a,e11.4)') shift//'ux1 ',bcdef%ux1
    write(*,'(a,e11.4)') shift//'vx0 ',bcdef%vx0
    write(*,'(a,e11.4)') shift//'vx1 ',bcdef%vx1
    write(*,'(a,e11.4)') shift//'wx0 ',bcdef%wx0
    write(*,'(a,e11.4)') shift//'wx1 ',bcdef%wx1
    write(*,'(a,e11.4)') shift//'uy0 ',bcdef%uy0
    write(*,'(a,e11.4)') shift//'uy1 ',bcdef%uy1
    write(*,'(a,e11.4)') shift//'vy0 ',bcdef%vy0
    write(*,'(a,e11.4)') shift//'vy1 ',bcdef%vy1
    write(*,'(a,e11.4)') shift//'wy0 ',bcdef%wy0
    write(*,'(a,e11.4)') shift//'wy1 ',bcdef%wy1
    write(*,'(a,e11.4)') shift//'uz0 ',bcdef%uz0
    write(*,'(a,e11.4)') shift//'uz1 ',bcdef%uz1
    write(*,'(a,e11.4)') shift//'vz0 ',bcdef%vz0
    write(*,'(a,e11.4)') shift//'vz1 ',bcdef%vz1
    write(*,'(a,e11.4)') shift//'wz0 ',bcdef%wz0
    write(*,'(a,e11.4)') shift//'wz1 ',bcdef%wz1
  case('segmented_s_line')
    write(*,'(a)') shift//'Boundary condition parameters'
    do i=1,12
      write(*,'(a,i2,a,e11.4)') shift//'bc_parameter ',i,' ',bcdef%bc_parameters(i)
    enddo
  end select
  write(*,'(a,e11.4)') shift//'utrans ',bcdef%utrans
  write(*,'(a,e11.4)') shift//'vtrans ',bcdef%vtrans
  write(*,'(a,l1)') shift//'isostasy ',params%isostasy
  write(*,'(a,l1)') shift//'flexure ',params%flexure
  write(*,'(a,l1)') shift//'isobc ',params%isobc 
  write(*,'(a,e11.4)') shift//'elastic_plate_thickness ',params%elastic_plate_thickness
  write(*,'(a,e11.4)') shift//'density_difference ',params%density_difference
  write(*,'(a,l1)') shift//'nest ',params%nest
  if (params%nest) then
    write(*,'(a)')  shift//'--- Nest parameters ---'
    write(*,'(a,a)')  shift//'lsoutfile ',trim(nest%lsoutfile)
    write(*,'(a,e11.4)') shift//'sselemx ',nest%sselemx
    write(*,'(a,e11.4)') shift//'sselemy ',nest%sselemy
    write(*,'(a,e11.4)') shift//'sselemz ',nest%sselemz
    write(*,'(a,e11.4)') shift//'xminls ',nest%xminls
    write(*,'(a,e11.4)') shift//'yminls ',nest%yminls
    write(*,'(a,e11.4)') shift//'zminls ',nest%zminls
  endif
  write(*,'(a,e11.4)') shift//'distance_exponent ',params%distance_exponent
endif
if (params%debug.gt.1) then
  write(threadinfo%Logunit,'(a)')         'Input file values: '
  write(threadinfo%Logunit,'(a32,i4)')    'nstep ',params%nstep
  write(threadinfo%Logunit,'(a32,i4)')    'material0 ',material0
  write(threadinfo%Logunit,'(a32,l1)')    'bulkvisc ',params%bulkvisc
  write(threadinfo%Logunit,'(a32,l1)')    'init_e2d ',params%init_e2d
  write(threadinfo%Logunit,'(a32,e11.4)') 'e2d0 ',params%e2d0
  write(threadinfo%Logunit,'(a32,e11.4)') 'minp ',params%minp
  write(threadinfo%Logunit,'(a32,e11.4)') 'pressure0 ',params%pressure0
  do i=1,params%nmat
    write(threadinfo%Logunit,'(a,i4,a)')    '--- Properties for material ',i,' ---'
    write(threadinfo%Logunit,'(a32,e11.4)') 'density ',mat(i)%density
    write(threadinfo%Logunit,'(a32,e11.4)') 'viscosity ',mat(i)%viscosity
    write(threadinfo%Logunit,'(a32,e11.4)') 'penalty ',mat(i)%penalty
    write(threadinfo%Logunit,'(a32,e11.4)') 'expon ',mat(i)%expon
    write(threadinfo%Logunit,'(a32,e11.4)') 'activationenergy ',mat(i)%activationenergy
    write(threadinfo%Logunit,'(a32,e11.4)') 'expansion ',mat(i)%expansion
    write(threadinfo%Logunit,'(a32,e11.4)') 'diffusivity ',mat(i)%diffusivity
    write(threadinfo%Logunit,'(a32,e11.4)') 'heat ',mat(i)%heat
    write(threadinfo%Logunit,'(a32,a8)')    'plasticity_type ',mat(i)%plasticity_type
    if (trim(mat(i)%plasticity_type).ne.'No') then
      write(threadinfo%Logunit,'(a,i4)') 'Plasticity parameters for material ',i
      do j=1,9
        write(threadinfo%Logunit,'(a32,i1,a2,e11.4)') 'plasticity parameter ',j,': ',mat(i)%plasticity_parameters(j)
      enddo
    endif
    write (threadinfo%Logunit,'(a,i4)') 'Material transitions for material ',i
    do j=1,6
      write(threadinfo%Logunit,'(a32,i1,a2,e11.4)') 'mattrans ',j,' ',mat(i)%mattrans(j)
    enddo
    do j=1,6
      write(threadinfo%Logunit,'(a32,i1,a2,i4)') 'transnum ',j,' ',mat(i)%transnum(j)
    enddo
  enddo
  write(threadinfo%Logunit,'(a32,e11.4)') 'viscositymin ',params%viscositymin
  write(threadinfo%Logunit,'(a32,e11.4)') 'viscositymax ',params%viscositymax
  write(threadinfo%Logunit,'(a32,e11.4)') 'vex ',params%vex
  write(threadinfo%Logunit,'(a32,i4)')    'leveluniform_oct ',params%leveluniform_oct
  write(threadinfo%Logunit,'(a32,i4)')    'levelmax_oct ',params%levelmax_oct
  write(threadinfo%Logunit,'(a32,i4)')    'matrule ',params%matrule
  write(threadinfo%Logunit,'(a32,i4)')    'levelcut ',params%levelcut
  write(threadinfo%Logunit,'(a32,i4)')    'levelapprox ',params%levelapprox
  write(threadinfo%Logunit,'(a32,l1)')    'calculate_temp ',params%calculate_temp
  write(threadinfo%Logunit,'(a32,e11.4)') 'ztemp ',params%ztemp
  write(threadinfo%Logunit,'(a32,i4)')    'smoothing_type ',params%smoothing_type
  write(threadinfo%Logunit,'(a32,l1)')    'normaladvect ',params%normaladvect
  write(threadinfo%Logunit,'(a32,l1)')    'excl_vol ',params%excl_vol
  write(threadinfo%Logunit,'(a32,l1)')    'adaptive_tol ',params%adaptive_tol
  do i=1,params%ns
    write(threadinfo%Logunit,'(a,i4,a)')         '--- Properties for surface ',i,' ---'
    write(threadinfo%Logunit,'(a32,i4,a,i4)')    'surface ',i,' itype ',surface(i)%itype
    write(threadinfo%Logunit,'(a32,i4,a,i4)')    'surface ',i,' material ',surface(i)%material
    write(threadinfo%Logunit,'(a32,i4,a,l1)')    'surface ',i,' rand ',surface(i)%rand
    write(threadinfo%Logunit,'(a32,i4,a,i4)')    'surface ',i,' levelt ',surface(i)%levelt
    write(threadinfo%Logunit,'(a32,i4,a,e11.4)') 'surface ',i,' stretch ',surface(i)%stretch
    write(threadinfo%Logunit,'(a32,i4,a,i4)')    'surface ',i,' criterion ',surface(i)%criterion
    write(threadinfo%Logunit,'(a32,i4,a,e11.4)') 'surface ',i,' anglemax ',surface(i)%anglemax
    write(threadinfo%Logunit,'(a32,i4,a,e11.4)') 'surface ',i,' anglemaxoctree ',surface(i)%anglemaxoctree
    write(threadinfo%Logunit,'(a32,i4,a,i4)')    'surface ',i,' spread_surface_points ',surface(i)%spread_surface_points
    write(threadinfo%Logunit,'(a32,i4,a,i4)')    'surface ',i,' type ',surface(i)%surface_type
    write(threadinfo%Logunit,'(a32,i4,a,e11.4)') 'surface ',i,' sp01 ',surface(i)%sp01
    write(threadinfo%Logunit,'(a32,i4,a,e11.4)') 'surface ',i,' sp02 ',surface(i)%sp02
    write(threadinfo%Logunit,'(a32,i4,a,e11.4)') 'surface ',i,' sp03 ',surface(i)%sp03
    write(threadinfo%Logunit,'(a32,i4,a,e11.4)') 'surface ',i,' sp04 ',surface(i)%sp04
    write(threadinfo%Logunit,'(a32,i4,a,e11.4)') 'surface ',i,' sp05 ',surface(i)%sp05
    write(threadinfo%Logunit,'(a32,i4,a,e11.4)') 'surface ',i,' sp06 ',surface(i)%sp06
    write(threadinfo%Logunit,'(a32,i4,a,e11.4)') 'surface ',i,' sp07 ',surface(i)%sp07
    write(threadinfo%Logunit,'(a32,i4,a,e11.4)') 'surface ',i,' sp08 ',surface(i)%sp08
    write(threadinfo%Logunit,'(a32,i4,a,e11.4)') 'surface ',i,' sp09 ',surface(i)%sp09
    write(threadinfo%Logunit,'(a32,i4,a,e11.4)') 'surface ',i,' sp10 ',surface(i)%sp10
    write(threadinfo%Logunit,'(a32,i4,a,e11.4)') 'surface ',i,' sp11 ',surface(i)%sp11
    write(threadinfo%Logunit,'(a32,i4,a,e11.4)') 'surface ',i,' sp12 ',surface(i)%sp12
    write(threadinfo%Logunit,'(a32,i4,a,e11.4)') 'surface ',i,' activation_time ',surface(i)%activation_time
    write(threadinfo%Logunit,'(a32,i4,a,i4)')    'surface ',i,' leveloct ',surface(i)%leveloct
  enddo
  write(threadinfo%Logunit,'(a32,l1)') 'ref_on_faces ',params%ref_on_faces
  if (params%ref_on_faces) then
    do i=1,6
      write(threadinfo%Logunit,'(a,i1,a)')    '--- Refinement values for face ',i,' ---'
      write(threadinfo%Logunit,'(a32,i4)')    'cube_faces level ',cube_faces(i)%level
      write(threadinfo%Logunit,'(a32,e11.4)') 'cube_faces l ',cube_faces(i)%l
      write(threadinfo%Logunit,'(a32,e11.4)') 'cube_faces r ',cube_faces(i)%r
      write(threadinfo%Logunit,'(a32,e11.4)') 'cube_faces b ',cube_faces(i)%b
      write(threadinfo%Logunit,'(a32,e11.4)') 'cube_faces t ',cube_faces(i)%t
    enddo
  endif
  write(threadinfo%Logunit,'(a32,i12)') 'noctreemax ',params%noctreemax
  write(threadinfo%Logunit,'(a32,i4)') 'nonlinear_iterations ',params%nonlinear_iterations
  write(threadinfo%Logunit,'(a32,i4)') 'initial_refine_level ',params%initial_refine_level
  write(threadinfo%Logunit,'(a32,e11.4)') 'dt ',params%dt
  write(threadinfo%Logunit,'(a32,l1)') 'damp_surface ',params%damp_surface
  write(threadinfo%Logunit,'(a32,e11.4)') 'damp_factor ',params%damp_factor
  write(threadinfo%Logunit,'(a32,l1)') 'remove_surf_pts ',params%remove_surf_pts
  write(threadinfo%Logunit,'(a32,e11.4)') 'penalty ',params%penalty
  write(threadinfo%Logunit,'(a32,e11.4)') 'tempscale ',params%tempscale
  write(threadinfo%Logunit,'(a32,e11.4)') 'refine_ratio ',params%refine_ratio
  write(threadinfo%Logunit,'(a32,i4)') 'refine_criterion ',params%refine_criterion
  write(threadinfo%Logunit,'(a32,e11.4)') 'octree_refine_ratio ',params%octree_refine_ratio
  write(threadinfo%Logunit,'(a32,e11.4)') 'courant ',params%courant
  write(threadinfo%Logunit,'(a32,l1)') 'invariants_2d ',params%invariants_2d
  write(threadinfo%Logunit,'(a32,i4)') 'npmin ',params%npmin
  write(threadinfo%Logunit,'(a32,i4)') 'npmax ',params%npmax
  write(threadinfo%Logunit,'(a32,i4)') 'griditer ',params%griditer
  write(threadinfo%Logunit,'(a32,e11.4)') 'tol ',params%tol
  write(threadinfo%Logunit,'(a32,i4)') 'niter_move ',params%niter_move
  write(threadinfo%Logunit,'(a32,l1)') 'ismooth ',params%ismooth
  write(threadinfo%Logunit,'(a32,i4)') 'nb_iter_nl_min ',params%nb_iter_nl_min
  write(threadinfo%Logunit,'(a32,l1)') 'visualise_matrix ',params%visualise_matrix
  write(threadinfo%Logunit,'(a32,l1)') 'renumber_nodes ',params%renumber_nodes
  if (params%nboxes.gt.0) then
    do i=1,params%nboxes
      write(threadinfo%Logunit,'(a,i3,a)')    '--- Input values for box ',i,' ---'
      write(threadinfo%Logunit,'(a32,e11.4)') 'box x0 ',boxes(i)%x0
      write(threadinfo%Logunit,'(a32,e11.4)') 'box x1 ',boxes(i)%x1
      write(threadinfo%Logunit,'(a32,e11.4)') 'box y0 ',boxes(i)%y0
      write(threadinfo%Logunit,'(a32,e11.4)') 'box y1 ',boxes(i)%y1
      write(threadinfo%Logunit,'(a32,e11.4)') 'box z0 ',boxes(i)%z0
      write(threadinfo%Logunit,'(a32,e11.4)') 'box z1 ',boxes(i)%z1
      write(threadinfo%Logunit,'(a32,i4)')    'box level ',boxes(i)%level
    enddo
  endif
  if (params%nsections.gt.0) then
    do i=1,params%nsections
      write(threadinfo%Logunit,'(a,i3,a)')    '--- Input values for section ',i,' ---'    
      write(threadinfo%Logunit,'(a32,i4)')    'section xyz ',sections(i)%xyz
      write(threadinfo%Logunit,'(a32,e11.4)') 'section slice ',sections(i)%slice
      write(threadinfo%Logunit,'(a32,l1)')    'section flag press ',sections(i)%flag_press
      write(threadinfo%Logunit,'(a32,l1)')    'section flag spress ',sections(i)%flag_spress
      write(threadinfo%Logunit,'(a32,l1)')    'section flag e2d ',sections(i)%flag_e2d
      write(threadinfo%Logunit,'(a32,l1)')    'section flag e3d ',sections(i)%flag_e3d
      write(threadinfo%Logunit,'(a32,l1)')    'section flag strain ',sections(i)%flag_strain
      write(threadinfo%Logunit,'(a32,l1)')    'section flag crit ',sections(i)%flag_crit
      write(threadinfo%Logunit,'(a32,l1)')    'section flag grid ',sections(i)%flag_grid
      write(threadinfo%Logunit,'(a32,l1)')    'section flag u ',sections(i)%flag_u
      write(threadinfo%Logunit,'(a32,l1)')    'section flag v ',sections(i)%flag_v  
      write(threadinfo%Logunit,'(a32,l1)')    'section flag w ',sections(i)%flag_w 
      write(threadinfo%Logunit,'(a32,l1)')    'section flag uvw ',sections(i)%flag_uvw
      write(threadinfo%Logunit,'(a32,l1)')    'section flag colour ',sections(i)%flag_colour
      write(threadinfo%Logunit,'(a32,l1)')    'section flag mu ',sections(i)%flag_mu
      write(threadinfo%Logunit,'(a32,l1)')    'section flag plastic ',sections(i)%flag_plastic
      write(threadinfo%Logunit,'(a32,l1)')    'section flag q ',sections(i)%flag_q
      write(threadinfo%Logunit,'(a32,l1)')    'section flag lode ',sections(i)%flag_lode
      write(threadinfo%Logunit,'(a32,l1)')    'section flag vfield ',sections(i)%flag_vfield
      write(threadinfo%Logunit,'(a32,l1)')    'section flag lsf ',sections(i)%flag_lsf
      write(threadinfo%Logunit,'(a32,i3)')    'section flag temp ',sections(i)%flag_temp
      write(threadinfo%Logunit,'(a32,l1)')    'section flag velvect ',sections(i)%flag_velvect
      write(threadinfo%Logunit,'(a32,e11.4)') 'section scale ',sections(i)%scale
      write(threadinfo%Logunit,'(a32,a3)')    'section colormap ',sections(i)%colormap
      write(threadinfo%Logunit,'(a32,i4)')    'section ncolours ',sections(i)%ncolours
    enddo
  endif
  write(threadinfo%Logunit,'(a32,l1)') 'erosion ',params%erosion
  write(threadinfo%Logunit,'(a32,e11.4)') 'zerosion ',params%zerosion
  write(threadinfo%Logunit,'(a32,e11.4)') 'length_scale ',params%length_scale
  write(threadinfo%Logunit,'(a32,e11.4)') 'velocity_scale ',params%velocity_scale
  write(threadinfo%Logunit,'(a32,e11.4)') 'fluvial_erosion ',params%fluvial_erosion
  write(threadinfo%Logunit,'(a32,e11.4)') 'diffusion_erosion ',params%diffusion_erosion
  write(threadinfo%Logunit,'(a32,i4)') 'baselevelx0 ',params%baselevelx0
  write(threadinfo%Logunit,'(a32,i4)') 'baselevelx1 ',params%baselevelx1
  write(threadinfo%Logunit,'(a32,i4)') 'baselevely0 ',params%baselevely0
  write(threadinfo%Logunit,'(a32,i4)') 'baselevely1 ',params%baselevely1
  write(threadinfo%Logunit,'(a32,l1)') 'compute_qpgram ',params%compute_qpgram
  write(threadinfo%Logunit,'(a)') '--- Boundary condition parameters ---'
  write(threadinfo%Logunit,'(a32,l1)') 'invariants_2d ',params%invariants_2d
  write(threadinfo%Logunit,'(a32,a)') 'bctype ',bcdef%bctype
  select case(trim(bcdef%bctype))
  case('basic')
    write(threadinfo%Logunit,'(a32,a)') 'bcorder ',bcdef%bcorder
    write(threadinfo%Logunit,'(a32,l1)') 'fixux0 ',bcdef%fixux0
    write(threadinfo%Logunit,'(a32,l1)') 'fixux1 ',bcdef%fixux1
    write(threadinfo%Logunit,'(a32,l1)') 'fixvx0 ',bcdef%fixvx0
    write(threadinfo%Logunit,'(a32,l1)') 'fixvx1 ',bcdef%fixvx1
    write(threadinfo%Logunit,'(a32,l1)') 'fixwx0 ',bcdef%fixwx0
    write(threadinfo%Logunit,'(a32,l1)') 'fixwx1 ',bcdef%fixwx1
    write(threadinfo%Logunit,'(a32,l1)') 'fixuy0 ',bcdef%fixuy0
    write(threadinfo%Logunit,'(a32,l1)') 'fixuy1 ',bcdef%fixuy1
    write(threadinfo%Logunit,'(a32,l1)') 'fixvy0 ',bcdef%fixvy0
    write(threadinfo%Logunit,'(a32,l1)') 'fixvy1 ',bcdef%fixvy1
    write(threadinfo%Logunit,'(a32,l1)') 'fixwy0 ',bcdef%fixwy0
    write(threadinfo%Logunit,'(a32,l1)') 'fixwy1 ',bcdef%fixwy1
    write(threadinfo%Logunit,'(a32,l1)') 'fixuz0 ',bcdef%fixuz0
    write(threadinfo%Logunit,'(a32,l1)') 'fixuz1 ',bcdef%fixuz1
    write(threadinfo%Logunit,'(a32,l1)') 'fixvz0 ',bcdef%fixvz0
    write(threadinfo%Logunit,'(a32,l1)') 'fixvz1 ',bcdef%fixvz1
    write(threadinfo%Logunit,'(a32,l1)') 'fixwz0 ',bcdef%fixwz0
    write(threadinfo%Logunit,'(a32,l1)') 'fixwz1 ',bcdef%fixwz1
    write(threadinfo%Logunit,'(a32,e11.4)') 'ux0 ',bcdef%ux0
    write(threadinfo%Logunit,'(a32,e11.4)') 'ux1 ',bcdef%ux1
    write(threadinfo%Logunit,'(a32,e11.4)') 'vx0 ',bcdef%vx0
    write(threadinfo%Logunit,'(a32,e11.4)') 'vx1 ',bcdef%vx1
    write(threadinfo%Logunit,'(a32,e11.4)') 'wx0 ',bcdef%wx0
    write(threadinfo%Logunit,'(a32,e11.4)') 'wx1 ',bcdef%wx1
    write(threadinfo%Logunit,'(a32,e11.4)') 'uy0 ',bcdef%uy0
    write(threadinfo%Logunit,'(a32,e11.4)') 'uy1 ',bcdef%uy1
    write(threadinfo%Logunit,'(a32,e11.4)') 'vy0 ',bcdef%vy0
    write(threadinfo%Logunit,'(a32,e11.4)') 'vy1 ',bcdef%vy1
    write(threadinfo%Logunit,'(a32,e11.4)') 'wy0 ',bcdef%wy0
    write(threadinfo%Logunit,'(a32,e11.4)') 'wy1 ',bcdef%wy1
    write(threadinfo%Logunit,'(a32,e11.4)') 'uz0 ',bcdef%uz0
    write(threadinfo%Logunit,'(a32,e11.4)') 'uz1 ',bcdef%uz1
    write(threadinfo%Logunit,'(a32,e11.4)') 'vz0 ',bcdef%vz0
    write(threadinfo%Logunit,'(a32,e11.4)') 'vz1 ',bcdef%vz1
    write(threadinfo%Logunit,'(a32,e11.4)') 'wz0 ',bcdef%wz0
    write(threadinfo%Logunit,'(a32,e11.4)') 'wz1 ',bcdef%wz1
  case('segmented_s_line')
    write(threadinfo%Logunit,'(a)') 'Boundary condition parameters'
    do i=1,12
      write(threadinfo%Logunit,'(a32,i2,a,e11.4)') 'bc_parameter ',i,' ',bcdef%bc_parameters(i)
    enddo
  end select
  write(threadinfo%Logunit,'(a32,e11.4)') 'utrans ',bcdef%utrans
  write(threadinfo%Logunit,'(a32,e11.4)') 'vtrans ',bcdef%vtrans
  write(threadinfo%Logunit,'(a32,l1)') 'isostasy ',params%isostasy
  write(threadinfo%Logunit,'(a32,l1)') 'flexure ',params%flexure
  write(threadinfo%Logunit,'(a32,l1)') 'isobc ',params%isobc 
  write(threadinfo%Logunit,'(a32,e11.4)') 'elastic_plate_thickness ',params%elastic_plate_thickness
  write(threadinfo%Logunit,'(a32,e11.4)') 'density_difference ',params%density_difference
  write(threadinfo%Logunit,'(a32,l1)') 'nest ',params%nest
  if (params%nest) then
    write(threadinfo%Logunit,'(a)') '--- Nest parameters ---'
    write(threadinfo%Logunit,'(a32,a)') 'lsoutfile ',trim(nest%lsoutfile)
    write(threadinfo%Logunit,'(a32,e11.4)') 'sselemx ',nest%sselemx
    write(threadinfo%Logunit,'(a32,e11.4)') 'sselemy ',nest%sselemy
    write(threadinfo%Logunit,'(a32,e11.4)') 'sselemz ',nest%sselemz
    write(threadinfo%Logunit,'(a32,e11.4)') 'xminls ',nest%xminls
    write(threadinfo%Logunit,'(a32,e11.4)') 'yminls ',nest%yminls
    write(threadinfo%Logunit,'(a32,e11.4)') 'zminls ',nest%zminls
  endif
  write(threadinfo%Logunit,'(a32,e11.4)') 'distance_exponent ',params%distance_exponent
endif

if(iproc.eq.0) call flush(8)
if(params%debug.gt.1) call flush(threadinfo%logunit)

if (iproc.eq.0) call system ('rm fort.8')

end subroutine read_input_file
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|