define_bc_rot_sub_init.f90

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              DEFINE_BC_SEGMENTED_S_LINE   Feb. 2009                          |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

subroutine define_bc_rot_sub_init(params,osolve,vo,bcdef)

!------------------------------------------------------------------------------|
!((((((((((((((((  Purpose of the routine  )))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! This routine assigns the velocity boundary conditions for the rotational     |
! subduction geometry                                                          |
!------------------------------------------------------------------------------|
!((((((((((((((((  declaration of the subroutine arguments  ))))))))))))))))))))
!------------------------------------------------------------------------------|

use definitions
!use mpi

implicit none

include 'mpif.h' 

type (parameters) params
type (octreesolve) osolve
type (void) vo
type (bc_definition) bcdef
type (nest_info) nest

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

integer i,iproc,nproc,ierr,tempcase
double precision :: eps,pi,y0,h,rx,wy,vin
double precision :: rmain,zmain,omega,theta,distmain,distrot
double precision :: nl_heat,hr,depth,ztemp,shellwidth

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

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

!basic constants
eps=1.d-10
pi=atan(1.d0)*4.d0

!set fixed to 0 everywhere
osolve%kfix=0
osolve%kfixt=0


!parameters from input file 
y0=bcdef%bc_parameters(1)
h=bcdef%bc_parameters(2)
rx=bcdef%bc_parameters(3)
wy=bcdef%bc_parameters(4)
vin=bcdef%bc_parameters(5)
shellwidth = bcdef%bc_parameters(6)    
tempcase=idint(bcdef%bc_parameters(7))
nl_heat=bcdef%bc_parameters(8)/params%tempscale
hr=bcdef%bc_parameters(9)

ztemp = params%ztemp*params%vex

!Derived parameters from input values
rmain=(h**2+rx**2)/(2.d0*h)
zmain=h-rmain
omega=vin/rmain

!rotating element is sphere

do i=1,osolve%nnode

!distance from center of sphere
distmain=sqrt(osolve%x(i)**2+((osolve%y(i)-y0)*rx/(0.5d0*wy))**2 + &
              (osolve%z(i)*params%vex-zmain)**2)
!distance from axis of rotation
distrot=sqrt(osolve%x(i)**2+(osolve%z(i)*params%vex-zmain)**2)

!setting velocity in indenter
if (distmain<rmain) then
    theta=asin(osolve%x(i)/distrot)
    osolve%kfix((i-1)*3+1)=1
    osolve%u(i)=distrot*omega*cos(theta)
    osolve%kfix((i-1)*3+2)=1
    osolve%v(i)=0.d0
    osolve%kfix((i-1)*3+3)=1
    osolve%w(i)=-osolve%x(i)*omega
elseif (distmain > rmain+shellwidth) 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
end if

!y=0 face free slip 
if (osolve%y(i)<eps) then
    !no flux through face
    osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
end if

!y=1 face free slip
if (osolve%y(i)>(1.d0-eps)) then
    !no flux through face
    osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
end if

!x=0 face
if (osolve%x(i)<eps) then
    !set temperature
    osolve%kfixt(i) = 1

    depth = ztemp - osolve%z(i)*params%vex
    select case (tempcase)
        case (1)
        !"typical isotherm taken from R.Bendicks simulation
        !nl_heat is scaled nonlinear heat component hr^2*S_o/k_T
            osolve%temp(i) = (depth/ztemp) + & 
                nl_heat*((1.d0-exp(-depth/hr))-(1.d0-exp(-ztemp/hr))*depth/ztemp)
        case default
            osolve%temp(i) = depth/ztemp
    end select
end if

!z=0 face
if (osolve%z(i) < eps) then
    if (distmain > rmain+shellwidth) then
        osolve%kfixt(i) = 1
        osolve%temp(i) = 1.d0
    end if      
end if

!x=1 face with influx vback
!if (osolve%x(i)>1.d0-eps) then
!end if

!z=1 face 
if (osolve%z(i) > 1.d0-eps) then
    osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=0.d0 !x-direction
    osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0 !y-direction
    osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0 !z-direction
    osolve%kfixt(i)=1 ; osolve%temp(i)=0.d0
end if

if (.not.vo%influid(i)) then
   osolve%kfixt(i)=1
   osolve%temp(i)=0.d0
endif

end do

end

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