diff --git a/define_bc_DT065.f90 b/define_bc_DT065.f90
new file mode 100644
index 0000000000000000000000000000000000000000..94c09a2e1d77e149e81ceea1363e827c49f00095
--- /dev/null
+++ b/define_bc_DT065.f90
@@ -0,0 +1,191 @@
+!------------------------------------------------------------------------------|
+!------------------------------------------------------------------------------|
+! |
+! ||===\\ |
+! || \\ |
+! || || //==\\ || || //==|| ||/==\\ |
+! || || || || || || || || || || |
+! || // || || || || || || || |
+! ||===// \\==// \\==\\ \\==\\ || |
+! |
+!------------------------------------------------------------------------------|
+!------------------------------------------------------------------------------|
+! |
+! DEFINE_BC_TIBET Feb. 2009 |
+! |
+!------------------------------------------------------------------------------|
+!------------------------------------------------------------------------------|
+
+subroutine define_bc (params,osolve,vo)
+
+!------------------------------------------------------------------------------|
+!(((((((((((((((( 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
+
+!------------------------------------------------------------------------------|
+!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
+!------------------------------------------------------------------------------|
+
+integer i,iproc,nproc,ierr
+double precision eps,lsf0,pi,l,h,x0,x1,x2,phi,spu,spv
+integer ie,ij,j,jp
+
+!------------------------------------------------------------------------------|
+!------------------------------------------------------------------------------|
+
+
+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.2945d0
+h=0.0625d0
+x0=0.5d0
+phi=45.d0
+phi=phi/180.d0*pi
+x1=x0-l*(1.d0-cos(phi))/sin(phi)
+x2=x0+l/tan(phi)
+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')
+
+ do i=1,osolve%nnode
+ 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.0.575d0) osolve%u(i)=(1.d0-sin(pi*(osolve%y(i)+.425d0)/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.0.575d0) osolve%u(i)=cos(phi)*(1.d0-sin(pi*(osolve%y(i)+.425d0)/2.d0)**100)
+ if (osolve%y(i).le.0.575d0) osolve%w(i)=-sin(phi)*(1.d0-sin(pi*(osolve%y(i)+.425d0)/2.d0)**100)
+ elseif (osolve%x(i).lt.x1) then
+ if (osolve%y(i).le.0.575d0) osolve%u(i)=(1.d0-sin(pi*(osolve%y(i)+.425d0)/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
+
+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
+!------------------------------------------------------------------------------|
+!------------------------------------------------------------------------------|