diff --git a/src/define_bc.f90 b/src/define_bc.f90
deleted file mode 100644
index b559cdfe31d965a535d49ff74b026be409c344d8..0000000000000000000000000000000000000000
--- a/src/define_bc.f90
+++ /dev/null
@@ -1,229 +0,0 @@
-!------------------------------------------------------------------------------|
-!------------------------------------------------------------------------------|
-! |
-! ||===\\ |
-! || \\ |
-! || || //==\\ || || //==|| ||/==\\ |
-! || || || || || || || || || || |
-! || // || || || || || || || |
-! ||===// \\==// \\==\\ \\==\\ || |
-! |
-!------------------------------------------------------------------------------|
-!------------------------------------------------------------------------------|
-! |
-! DEFINE_BC_TIBET Feb. 2009 |
-! |
-!------------------------------------------------------------------------------|
-!------------------------------------------------------------------------------|
-
-subroutine define_bc (params,osolve,vo,zisodisp)
-
-!------------------------------------------------------------------------------|
-!(((((((((((((((( Purpose of the routine )))))))))))))))))))))))))))))))))))))
-!------------------------------------------------------------------------------|
-! this routine assigns the boundary condition for the Tibet experiment
-
-!------------------------------------------------------------------------------|
-!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
-!------------------------------------------------------------------------------|
-
-use definitions
-
-implicit none
-
-type (parameters) params
-type (octreesolve) osolve
-type (void) vo
-double precision zisodisp(2**params%levelmax_oct+1,2**params%levelmax_oct+1)
-
-!------------------------------------------------------------------------------|
-!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
-!------------------------------------------------------------------------------|
-
-integer i,iproc,nproc,ierr
-double precision eps,lsf0,pi,lorig,h,x1,x2,phi,spu,spv,yend,cper,cscl,xoffset
-double precision ystart,theta,xmax,ykink,l,vin
-double precision,dimension(:),allocatable :: x0,ldisp
-integer ie,ij,j,jp
-logical firstcall
-
-!------------------------------------------------------------------------------|
-!------------------------------------------------------------------------------|
-
-INCLUDE 'mpif.h'
-
-call mpi_comm_size (mpi_comm_world,nproc,ierr)
-call mpi_comm_rank (mpi_comm_world,iproc,ierr)
-
-eps=1.d-10
-osolve%kfix=0
-osolve%kfixt=0
-pi=atan(1.d0)*4.d0
-l=0.15566d0 ! Mantle thickness (initial depth of s-line)
-ystart=0.2d0 ! Length of straight indenter interface
-ykink=0.3865d0 ! Position at which to start subduction zone kink
-yend=0.5d0 ! End of indenter
-theta=65.d0 ! Subduction zone kink angle
-theta=theta/180.d0*pi
-phi=45.d0 ! Subduction angle
-phi=phi/180.d0*pi
-xoffset=0.3d0 ! S-line curve y-intercept
-xmax=0.7d0
-vin=1.d0 ! Reference velocity
-spu=0.d0 ! x-velocity of s-point relative to incoming plate
-spv=0.d0 ! y-velocity of s-point relative to incoming plate
-
-select case(trim(params%infile))
-
-case('input.txt','input.small.txt')
- allocate(x0(osolve%nnode))
- do i=1,osolve%nnode
- if (osolve%y(i).le.ystart) then
- x0(i)=xoffset
- else if (osolve%y(i).le.ykink) then
- x0(i)=(osolve%y(i)-ystart)*tan(theta)+xoffset
- else
- x0(i)=xmax
- endif
- enddo
-
- if (params%isobc) then
- allocate (ldisp(osolve%nnode))
- firstcall=.true.
- call define_isostasy_bc(params,osolve,vo,zisodisp,firstcall,l,x0,vin,spu,& ! First call, to get depth of s-line
- spv,ldisp)
- endif
-
- do i=1,osolve%nnode
- if (params%isobc) l=ldisp(i)
- !x1=x0(i)-l*(1.d0-cos(phi))/sin(phi)
- x1=x0(i)-l/tan(phi)
- x2=x0(i)+l/tan(phi)
-
- if (osolve%x(i).lt.eps) then
- osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=0.d0
- osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
- osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
- if (osolve%y(i).le.(yend+0.075d0)) osolve%u(i)=(1.d0-sin(pi*(osolve%y(i)+(yend-0.075d0))/2.d0)**100)
- endif
- if (osolve%x(i).gt.1.d0-eps) then
- osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=0.d0
-! osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
-! osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
- endif
- if (osolve%y(i).lt.eps) then
-! osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=0.d0
- osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
-! osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
- endif
- if (osolve%y(i).gt.1.d0-eps) then
-! osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=0.d0
- osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
-! osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
- endif
- if (osolve%z(i).lt.eps) then
- osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=0.d0
- osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
- osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
- if ((osolve%x(i)-x1)*(osolve%x(i)-x2).le.0.d0) then
- if (osolve%y(i).le.(yend+0.075d0)) osolve%u(i)=cos(phi)*(1.d0-sin(pi*(osolve%y(i)+(yend-0.075d0))/2.d0)**100)
- if (osolve%y(i).le.(yend+0.075d0)) osolve%w(i)=-sin(phi)*(1.d0-sin(pi*(osolve%y(i)+(yend-0.075d0))/2.d0)**100)
- elseif (osolve%x(i).lt.x1) then
- if (osolve%y(i).le.(yend+0.075d0)) osolve%u(i)=(1.d0-sin(pi*(osolve%y(i)+(yend-0.075d0))/2.d0)**100)
- endif
- osolve%kfixt(i)=1
- osolve%temp(i)=1.d0
- endif
- if (osolve%z(i).gt.1.d0-eps) then
- osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=0.d0
- osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
- osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
- osolve%kfixt(i)=1
- osolve%temp(i)=0.d0
- endif
- if (.not.vo%influid(i)) then
- osolve%kfixt(i)=1
- osolve%temp(i)=0.d0
- endif
- if (spu.ne.0.d0) then ! If x-comp of s-point velocity is non-zero then
- osolve%u(i)=osolve%u(i)-spu ! subract s-point velocity from entire model domain
- endif ! to keep that point stationary
- if (spv.ne.0.d0) then ! Same condition as above in y-direction
- osolve%v(i)=osolve%v(i)-spv
- endif
- enddo
-
- if (params%isobc) then
- firstcall=.false.
- call define_isostasy_bc(params,osolve,vo,zisodisp,firstcall,l,x0,vin,spu,& ! Second call, to modify velocity B/Cs
- spv,ldisp)
- endif
- deallocate(x0)
- if(params%isobc) deallocate(ldisp)
-
-case ('input.jgr')
- call define_bc_jgr (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo)
-
-case ('input.sphere')
- call define_bc_sphere (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo)
-
-case ('input.3Dpunch')
- call define_bc_3Dpunch (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo)
-
-case ('input.2Dpunch')
- call define_bc_2Dpunch (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo)
-
-case ('input.2Dpunch_vert')
- call define_bc_2Dpunch_vert (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo)
-
-case ('input.folding')
- call define_bc_folding (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo)
-
-case ('input.pipo')
- call define_bc_pipo (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo)
-
-case ('input.ritske')
- call define_bc_ritske (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo)
-
-case ('input.ritske_isurf')
- call define_bc_ritske_isurf (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo)
-
-case ('input.subduction')
- call define_bc_subduction (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo)
-
-case ('input.riedel')
- call define_bc_riedel (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo)
-
-case ('input.parallipipede')
- call define_bc_parallipipede (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo)
-
-case default
- if (iproc.eq.0) print *,params%infile
- call stop_run (' pb with input file in define_bc.f90')
-
-end select
-
-! here we need to impose that the bc satisfy the bad faces constrains
-
-do ie=1,osolve%nface
- do j=1,4
- jp=1+mod(j,4)
- i=osolve%iface(j+4,ie)
- if (osolve%kfix((i-1)*3+1).eq.1) osolve%u(i)=(osolve%u(osolve%iface(j,ie))+osolve%u(osolve%iface(jp,ie)))/2.d0
- if (osolve%kfix((i-1)*3+2).eq.1) osolve%v(i)=(osolve%v(osolve%iface(j,ie))+osolve%v(osolve%iface(jp,ie)))/2.d0
- if (osolve%kfix((i-1)*3+3).eq.1) osolve%w(i)=(osolve%w(osolve%iface(j,ie))+osolve%w(osolve%iface(jp,ie)))/2.d0
- if (osolve%kfixt(i).eq.1) osolve%temp(i)=(osolve%temp(osolve%iface(j,ie))+osolve%temp(osolve%iface(jp,ie)))/2.d0
- enddo
-i=osolve%iface(9,ie)
-if (osolve%kfix((i-1)*3+1).eq.1) osolve%u(i)=sum(osolve%u(osolve%iface(1:4,ie)))/4.d0
-if (osolve%kfix((i-1)*3+2).eq.1) osolve%v(i)=sum(osolve%v(osolve%iface(1:4,ie)))/4.d0
-if (osolve%kfix((i-1)*3+3).eq.1) osolve%w(i)=sum(osolve%w(osolve%iface(1:4,ie)))/4.d0
-if (osolve%kfixt(i).eq.1) osolve%temp(i)=sum(osolve%temp(osolve%iface(1:4,ie)))/4.d0
-enddo
-
-if (params%debug==2) call output_bc (osolve)
-
-return
-end
-!------------------------------------------------------------------------------|
-!------------------------------------------------------------------------------|