!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! DEFINE_BC_SEGMENTED_S_LINE_ROUND Jan. 2014 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine define_bc_segmented_s_line_round(params,osolve,vo,bcdef)
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine )))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! This routine assigns the velocity boundary conditions for the rounded
! segmented s-line 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
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
integer i,iproc,nproc,ierr
double precision :: eps,pi,x1,x2,phi,yend,xstart,ystart
double precision :: theta,l,vin,vzfluxscl,cntvel,dxy,xend,xsym,ymax,xwidth
double precision :: ywidth,nb,uextra,vextra,arcradius,ystartarc
double precision :: yendarc
double precision,dimension(:),allocatable :: x0,ldisp
integer nelemx,nelemz,xdisp,ydisp
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
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=bcdef%bc_parameters(1)
phi=bcdef%bc_parameters(2)
theta=bcdef%bc_parameters(3)
xsym=bcdef%bc_parameters(4)
ystart=bcdef%bc_parameters(5)
yend=bcdef%bc_parameters(6)
ymax=bcdef%bc_parameters(7)
vin=bcdef%bc_parameters(8)
nelemx=idint(bcdef%bc_parameters(9))
nelemz=idint(bcdef%bc_parameters(10))
uextra=bcdef%bc_parameters(11)
vextra=bcdef%bc_parameters(12)
arcradius=bcdef%bc_parameters(13)
ywidth=yend-ystart
phi=phi/180.d0*pi
theta=theta/180.d0*pi
xwidth=ywidth*tan(theta)
xstart=xsym-(xwidth/2.d0)
xend=xsym+(xwidth/2.d0)
vzfluxscl=0.5d0*(1.d0+(tan(phi)/tan((pi-phi)/2.d0)))
nb=2**params%levelmax_oct
dxy=1.d0/2**(params%levelmax_oct+1.d0)
ystartarc=ystart-arcradius*tan(theta/2.d0) ! Define the y-start of the arcuate region
yendarc=ystart+arcradius*tan(theta/2.d0)*cos(theta) ! Define the y-end of the arcuate region
allocate(x0(osolve%nnode))
do i=1,osolve%nnode
if (osolve%y(i).le.ystartarc) then
x0(i)=xstart
else if (osolve%y(i).le.yendarc) then
x0(i)=xstart+arcradius-sqrt(arcradius**2.d0-(osolve%y(i)-ystartarc)**2.d0)
else if (osolve%y(i).le.yend) then
x0(i)=xstart+(osolve%y(i)-ystart)*tan(theta)
else
x0(i)=xend
endif
enddo
! write(*,*) 'arcradius: ',arcradius
! write(*,*) 'ystart: ',ystart
! write(*,*) 'ystartarc: ',ystartarc
! write(*,*) 'yendarc: ',yendarc
! write(*,*) 'x0: ',x0(:)
if (params%isobc) then
allocate (ldisp(osolve%nnode))
do i=1,osolve%nnode
xdisp=idnint(x0(i)*nb)+1
ydisp=idnint(osolve%y(i)*nb)+1
ldisp(i)=l+bcdef%zisodisp(xdisp,ydisp)
enddo
endif
do i=1,osolve%nnode
cntvel=0.d0
if (params%isobc) l=ldisp(i)
x1=x0(i)-l/tan(phi)
x2=x0(i)+l/tan(phi)
if (((x2-x0(i))-(x0(i)-x1)).gt.eps) vzfluxscl=1.d0 ! Set vzfluxscl equal to one if B/Cs are not symmetric
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.(ymax-(real(nelemx)*dxy))) then
osolve%u(i)=1.d0+uextra
elseif (osolve%y(i).le.(ymax+(real(nelemx)*dxy))) then
osolve%u(i)=(1.d0+uextra)*(1.d0-(osolve%y(i)-(ymax-(real(nelemx)*dxy)))/(2.d0*real(nelemx)*dxy))
endif
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).le.(x1-real(nelemz)*dxy)) then
if (osolve%y(i).le.(ymax-(real(nelemx)*dxy))) then
osolve%u(i)=1.d0+uextra
elseif (osolve%y(i).le.(ymax+(real(nelemx)*dxy))) then
osolve%u(i)=(1.d0+uextra)*(1.d0-(osolve%y(i)-(ymax-(real(nelemx)*dxy)))/(2.d0*real(nelemx)*dxy))
endif
elseif (osolve%x(i).le.(x1+real(nelemz)*dxy)) then
if (osolve%y(i).le.(ymax-(real(nelemx)*dxy))) then
osolve%u(i)=1.d0+uextra
elseif (osolve%y(i).le.(ymax+(real(nelemx)*dxy))) then
osolve%u(i)=(1.d0+uextra)*(1.d0-(osolve%y(i)-(ymax-(real(nelemx)*dxy)))/(2.d0*real(nelemx)*dxy))
endif
osolve%w(i)=(osolve%u(i)/(1.d0+uextra))*(osolve%x(i)-(x1-(real(nelemz)*dxy)))/(2.d0*real(nelemz)*dxy)*(vzfluxscl*(-sin(phi)))
osolve%u(i)=osolve%u(i)*(((1.d0-(osolve%x(i)-(x1-(real(nelemz)*dxy)))/(2.d0*real(nelemz)*dxy))*(1.d0-vzfluxscl*cos(phi)))+vzfluxscl*cos(phi))
elseif (osolve%x(i).le.(x2-real(nelemz)*dxy)) then
if (osolve%y(i).le.(ymax-(real(nelemx)*dxy))) then
osolve%u(i)=1.d0+uextra
elseif (osolve%y(i).le.(ymax+(real(nelemx)*dxy))) then
osolve%u(i)=(1.d0+uextra)*(1.d0-(osolve%y(i)-(ymax-(real(nelemx)*dxy)))/(2.d0*real(nelemx)*dxy))
endif
osolve%w(i)=(osolve%u(i)/(1.d0+uextra))*vzfluxscl*(-sin(phi))
osolve%u(i)=osolve%u(i)*vzfluxscl*cos(phi)
elseif (osolve%x(i).le.(x2+real(nelemz)*dxy)) then
if (osolve%y(i).le.(ymax-(real(nelemx)*dxy))) then
osolve%u(i)=1.d0+uextra
elseif (osolve%y(i).le.(ymax+(real(nelemx)*dxy))) then
osolve%u(i)=(1.d0+uextra)*(1.d0-(osolve%y(i)-(ymax-(real(nelemx)*dxy)))/(2.d0*real(nelemx)*dxy))
endif
osolve%w(i)=(osolve%u(i)/(1.d0+uextra))*((1.d0-(osolve%x(i)-(x2-(real(nelemz)*dxy)))/(2.d0*real(nelemz)*dxy))*(vzfluxscl*(-sin(phi))))
osolve%u(i)=osolve%u(i)*((1.d0-(osolve%x(i)-(x2-(real(nelemz)*dxy)))/(2.d0*real(nelemz)*dxy))*(vzfluxscl*cos(phi)))
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 (abs(bcdef%utrans).gt.eps) then
osolve%u(i)=osolve%u(i)+bcdef%utrans
endif
if (abs(bcdef%vtrans).gt.eps) then
osolve%v(i)=osolve%v(i)+bcdef%vtrans
endif
osolve%u(i)=osolve%u(i)*vin
osolve%v(i)=osolve%v(i)*vin
osolve%w(i)=osolve%w(i)*vin
enddo
if (params%isobc) then
call define_isostasy_bc(params,osolve,bcdef)
deallocate(ldisp)
endif
deallocate(x0)
end
!------------------------------------------------------------------------------|