!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! DEFINE_BC_SEGMENTED_S_LINE_PARABOLA August 2015 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine define_bc_segmented_s_line_parabola(params,osolve,vo,bcdef,nest)
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine )))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! This routine assigns the velocity boundary conditions for the 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
type (nest_info) nest
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
integer i,iproc,nproc,ierr
double precision :: eps,lsf0,pi,lorig,h,x1,x2,phi,yend,cper,cscl,xstart,ystart
double precision :: theta,l,vin,vzfluxscl,cntvel,dxy,xend,xsym,ymax,xwidth
double precision :: ywidth,xdisp,ydisp,nb
double precision,dimension(:),allocatable :: x0,ldisp
integer ie,ij,j,jp,nelemx,nelemz
double precision :: base,startp,endp,alpha,kink1,kink2,flytt,a,d
double precision :: wmax,uend,vmag,dipangle
double precision e(1),f(1),y(1)
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
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)
base=bcdef%bc_parameters(2)
startp=bcdef%bc_parameters(3)
endp=bcdef%bc_parameters(4)
vin=bcdef%bc_parameters(5)
alpha=bcdef%bc_parameters(6)
kink1=bcdef%bc_parameters(7)
kink2=bcdef%bc_parameters(8)
nelemx=idint(bcdef%bc_parameters(9))
nelemz=idint(bcdef%bc_parameters(10))
nb=2**params%levelmax_oct
dxy=1.d0/2**(params%levelmax_oct+1.d0)
alpha=alpha*pi/180
flytt=(kink2-kink1) * tan(alpha)
a = (base-l) /((endp-startp)*(endp-startp))
d = startp
wmax=(2*(vin+bcdef%utrans)*l)/(endp-startp)
uend=wmax/(2*a*(endp-startp))
vmag=sqrt(uend**2.d0+wmax**2.d0)
do i=1,osolve%nnode
if (osolve%x(i).lt.eps) then
osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=1.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%x(i).gt.1.d0-eps) then
osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=0.d0
endif
if (osolve%y(i).lt.eps) then
osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
endif
if (osolve%y(i).gt.1.d0-eps) then
osolve%kfix((i-1)*3+2)=1 ; osolve%v(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
osolve%kfixt(i)=1 ; osolve%temp(i)=1.d0
! Velocity in first margin-normal convergence segment
if (osolve%y(i) .lt. kink1) then
! Slab dip angle
dipangle = atan(2*a*(osolve%x(i)-startp))
! Velocity before transition to subduction
if (osolve%x(i) .le. startp-real(nelemz)*dxy) then
osolve%u(i) = vin
! Velocity in transition to subduction
elseif (osolve%x(i) .le. startp+real(nelemz)*dxy) then
! Velocity contribution from incoming horizontal velocity
! (diminishes to zero across transition)
osolve%u(i) = vin * (1.d0 - (osolve%x(i)-(startp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
! Velocity contribution from subduction region
! (increases from zero to full magnitude across transition)
osolve%u(i) = osolve%u(i) + vmag*cos(dipangle) * (osolve%x(i)-(startp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy)
osolve%w(i) = vmag*sin(dipangle) * (osolve%x(i)-(startp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy)
! Velocity in subduction region
elseif (osolve%x(i) .le. endp-real(nelemz)*dxy) then
osolve%u(i) = vmag*cos(dipangle)
osolve%w(i) = vmag*sin(dipangle)
! Velocity in transition from subduction region
elseif (osolve%x(i) .le. endp+real(nelemz)*dxy) then
! Velocity contribution from subduction region
! (decreases to zero across transition)
osolve%u(i) = vmag*cos(dipangle) * (1.d0 - (osolve%x(i)-(endp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
osolve%w(i) = vmag*sin(dipangle) * (1.d0 - (osolve%x(i)-(endp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
endif
! Velocity in oblique convergence segment
elseif (osolve%y(i) .lt. kink2) then
! Slab dip angle
dipangle = atan(2*a*(osolve%x(i)-((osolve%y(i)-kink1)*tan(alpha)+startp)))
! Velocity before transition to subduction
if (osolve%x(i) .le. (osolve%y(i)-kink1)*tan(alpha)+startp-real(nelemz)*dxy) then
osolve%u(i) = vin
! Velocity in transition to subduction
elseif (osolve%x(i) .le. (osolve%y(i)-kink1)*tan(alpha)+startp+real(nelemz)*dxy) then
! Velocity contribution from incoming horizontal velocity
! (diminishes to zero across transition)
osolve%u(i)= vin * (1.d0 - (osolve%x(i)-((osolve%y(i)-kink1)*tan(alpha)+startp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
! Velocity contribution from subduction region
! (increases from zero to full magnitude across transition)
osolve%u(i)= osolve%u(i) + vmag*cos(dipangle) * (osolve%x(i)-((osolve%y(i)-kink1)*tan(alpha)+startp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy)
osolve%w(i) = vmag*sin(dipangle) * (osolve%x(i)-((osolve%y(i)-kink1)*tan(alpha)+startp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy)
! Velocity in subduction region
elseif (osolve%x(i) .le. (osolve%y(i)-kink1)*tan(alpha)+endp-real(nelemz)*dxy) then
osolve%u(i) = vmag*cos(dipangle)
osolve%w(i) = vmag*sin(dipangle)
! Velocity in transition from subduction region
elseif (osolve%x(i) .le. (osolve%y(i)-kink1)*tan(alpha)+endp+real(nelemz)*dxy) then
! Velocity contribution from subduction region
! (decreases to zero across transition)
osolve%u(i) = vmag*cos(dipangle) * (1.d0 - (osolve%x(i)-((osolve%y(i)-kink1)*tan(alpha)+endp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
osolve%w(i) = vmag*sin(dipangle) * (1.d0 - (osolve%x(i)-((osolve%y(i)-kink1)*tan(alpha)+endp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
endif
! Velocity in second margin-normal convergence segment
else
! Slab dip angle
dipangle = atan(2*a*(osolve%x(i)-(flytt+startp)))
! Velocity before transition to subduction
if (osolve%x(i) .le. flytt+startp-real(nelemz)*dxy) then
osolve%u(i) = vin
! Velocity in transition to subduction
elseif (osolve%x(i) .le. flytt+startp+real(nelemz)*dxy) then
! Velocity contribution from incoming horizontal velocity
! (diminishes to zero across transition)
osolve%u(i) = vin * (1.d0 - (osolve%x(i)-((flytt+startp)-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
! Velocity contribution from subduction region
! (increases from zero to full magnitude across transition)
osolve%u(i) = osolve%u(i) + vmag * cos(dipangle) * (osolve%x(i)-((flytt+startp)-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy)
osolve%w(i) = vmag * sin(dipangle) * (osolve%x(i)-((flytt+startp)-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy)
! Velocity in subduction region
elseif (osolve%x(i) .le. flytt+endp-real(nelemz)*dxy) then
osolve%u(i) = vmag * cos(dipangle)
osolve%w(i) = vmag * sin(dipangle)
! Velocity in transition from subduction region
elseif (osolve%x(i) .le. flytt+endp+real(nelemz)*dxy) then
! Velocity contribution from subduction region
! (decreases to zero across transition)
osolve%u(i) = vmag * cos(dipangle) * (1.d0 - (osolve%x(i)-((flytt+endp)-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
osolve%w(i) = vmag * sin(dipangle) * (1.d0 - (osolve%x(i)-((flytt+endp)-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
endif
endif
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
if (osolve%kfix((i-1)*3+1)==1) osolve%u(i)=osolve%u(i)+bcdef%utrans
endif
if (abs(bcdef%vtrans).gt.eps) then
if (osolve%kfix((i-1)*3+2)==1) osolve%v(i)=osolve%v(i)+bcdef%vtrans
endif
if (osolve%kfix((i-1)*3+1)==1) osolve%u(i)=osolve%u(i)*vin
if (osolve%kfix((i-1)*3+2)==1) osolve%v(i)=osolve%v(i)*vin
if (osolve%kfix((i-1)*3+3)==1) osolve%w(i)=osolve%w(i)*vin
enddo
if (params%isobc) then
call define_isostasy_bc(params,osolve,bcdef)
endif
end
!------------------------------------------------------------------------------|