!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! DEFINE_BC_SEGMENTED_S_LINE Feb. 2009 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine define_bc_rot_subduction(params,osolve,vo,bcdef,nest)
!------------------------------------------------------------------------------|
!(((((((((((((((( 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,strictflow
double precision :: eps,pi,y0,h,rx,wy,vin,vtop,vback,nelemx,nelemz,dxy,dz
double precision :: rmain,zmain,omega,theta,phi,distmain,distrot,distz0plane
double precision :: distzonerad,distzoney,distzone,distedge
double precision :: xedge,xminedge,xmaxedge,yfix_xedge, yedge,yminedge,ymaxedge
double precision :: zedge,zminedge,zmaxedge, ustart,ustart2,uend,uend2
double precision :: wstart,wstart2,wend, temp2,tempstart,tempstart2,tempend
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
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)
vtop=bcdef%bc_parameters(6)
vback=-bcdef%bc_parameters(7) !change in orientation
nelemx=bcdef%bc_parameters(8)
nelemz=bcdef%bc_parameters(9)
tempcase=idint(bcdef%bc_parameters(10))
temp2=bcdef%bc_parameters(11)/params%tempscale
strictflow=idint(bcdef%bc_parameters(12))
dxy=nelemx/2.d0**(params%levelmax_oct+1.d0) !0.5*transition width in xy
dz=nelemz*params%vex/2.d0**(params%levelmax_oct+1.d0) !0.5*transition width in z
!Derived parameters from input values
rmain=(h**2+rx**2)/(2.d0*h)
zmain=h-rmain
omega=vin/rmain
!rotating element is sphere or cylinder
!case cylinder
if(wy<0) then
wy=-wy
!setting boundary conditions
do i=1,osolve%nnode
!distance from cylinder axis
distmain=sqrt(osolve%x(i)**2+(osolve%z(i)*params%vex-zmain)**2)
!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
if ((y0-0.5d0*wy) < eps) then
!rotating cylinder cuts y=0 face
if (distmain <= rmain) then
!angle of rotation
theta=asin(osolve%x(i)/distmain)
!set x velocity
osolve%kfix((i-1)*3+1)=1
osolve%u(i)=distmain*omega*cos(theta)
!set z velocity
osolve%kfix((i-1)*3+3)=1
osolve%w(i)=-osolve%x(i)*omega
end if
end if
end if
!y=1 face free slip
if (osolve%y(i)>(1.d0-eps)) then
osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
if ((y0+0.5d0*wy) > (1.d0-eps)) then
!rotating cylinder cuts y=1 face
if (distmain <= rmain) then
!angle of rotation
theta=asin(osolve%x(i)/distmain)
!set x velocity
osolve%kfix((i-1)*3+1)=1
osolve%u(i)=distmain*omega*cos(theta)
!set z velocity
osolve%kfix((i-1)*3+3)=1
osolve%w(i)=-osolve%x(i)*omega
end if
end if
end if
!x=0 face
if (osolve%x(i)<eps) then
!set y and z velocity
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
!transition zone in y-direction >1 for outside, <0 for inside
distzone=(abs(osolve%y(i)-y0)-0.5d0*wy+dxy)/(2.d0*dxy)
if (distzone > 1.d0) then
!outside rotating cylinder y-wise
osolve%kfix((i-1)*3+1)=1
osolve%u(i)=vback !x-direction
elseif (distzone < 0.d0) then
!inside rotating cylinder y-wise
if (osolve%z(i)*params%vex>(h+dz)) then
!above rotating cylinder z-wise
osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=vtop
elseif (osolve%z(i)*params%vex<(h-dz)) then
!within rotating cylinder z-wise
osolve%kfix((i-1)*3+1)=1
osolve%u(i)=omega*(osolve%z(i)*params%vex-zmain)
else
!transition zone z-wise
distzone=(osolve%z(i)*params%vex-(h-dz))/(2.d0*dz)
ustart=omega*(h-dz-zmain)
uend=vtop
osolve%kfix((i-1)*3+1)=1
osolve%u(i)=ustart+distzone*(uend-ustart)
end if
else
!transition zone y-wise at both sides (symmetry)
if (osolve%z(i)*params%vex>(h+dz)) then
!above rotating cylinder z-wise
osolve%kfix((i-1)*3+1)=1
ustart=vtop
uend=vback
osolve%u(i)=ustart+distzone*(uend-ustart)
elseif (osolve%z(i)*params%vex<(h-dz)) then
!within rotating cylinder z-wise
ustart=omega*(osolve%z(i)*params%vex-zmain)
uend=vback
osolve%kfix((i-1)*3+1)=1
osolve%u(i)=ustart+distzone*(uend-ustart)
else
!transition at corner
!start velocity scales according to z-transition zone
ustart2=omega*(h-dz-zmain)
ustart=ustart2+(osolve%z(i)*params%vex-(h-dz))/(2.d0*dz)*(vtop-ustart2)
!y transition to vback
uend=vback
osolve%kfix((i-1)*3+1)=1
osolve%u(i)=ustart+distzone*(uend-ustart)
end if
end if
end if
!z=0 face
if (osolve%z(i)*params%vex<eps) then
!set velocity
osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0 !y-direction
!transition zone in y-direction >1 for outside, <0 for inside
distzone=(abs(osolve%y(i)-y0)-0.5d0*wy+dxy)/(2.d0*dxy)
if (distzone > 1.d0) then
!outside rotating cylinder y-wise
osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=vback !x-direction
osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0 !z-direction
select case(tempcase)
case (0,2,3,4)
osolve%kfixt(i)=1
osolve%temp(i)=1.d0
case (1)
osolve%kfixt(i)=1
osolve%temp(i)=osolve%x(i)*(1.0d0-temp2)+temp2
case default
osolve%kfixt(i)=1
osolve%temp(i)=1.d0
end select
elseif (distzone < 0.d0) then
!inside rotating cylinder y-wise
if (osolve%x(i)>(rx+dxy)) then
!outside rotating cylinder x-wise
osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=vback
osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
select case(tempcase)
case (0,2,3,4)
osolve%kfixt(i)=1
osolve%temp(i)=1.d0
case (1)
osolve%kfixt(i)=1
osolve%temp(i)=osolve%x(i)*(1.0d0-temp2)+temp2
case default
osolve%kfixt(i)=1
osolve%temp(i)=1.d0
end select
elseif (osolve%x(i)<(rx-dxy)) then
!within rotating cylinder x-wise
distmain=sqrt(osolve%x(i)**2+zmain**2)
theta=asin(osolve%x(i)/distmain)
!set x velocity
osolve%kfix((i-1)*3+1)=1
osolve%u(i)=distmain*omega*cos(theta)
!set z velocity
osolve%kfix((i-1)*3+3)=1
osolve%w(i)=-osolve%x(i)*omega
select case(tempcase)
case (0)
osolve%kfixt(i)=1
osolve%temp(i)=1.d0
case (1)
osolve%kfixt(i)=1
osolve%temp(i)=osolve%x(i)*(1.0d0-temp2)+temp2
case (2)
osolve%kfixt(i)=1
osolve%temp(i)=temp2
case (3,4)
osolve%kfixt(i)=1
osolve%temp(i)=osolve%x(i)/rx*(temp2-1.0d0)+1.0d0
case default
osolve%kfixt(i)=1
osolve%temp(i)=1.d0
end select
else
!transition zone x-wise
!redefine transition zone in x-direction
distzone=(osolve%x(i)-(rx-dxy))/(2.d0*dxy)
!values for start of transition zone
distmain=sqrt((rx-dxy)**2+zmain**2)
theta=asin((rx-dxy)/distmain)
ustart=distmain*omega*cos(theta)
uend=vback
wstart=-(rx-dxy)*omega
wend=0.d0
!set velocities
osolve%kfix((i-1)*3+1)=1
osolve%u(i)=ustart+distzone*(uend-ustart)
osolve%kfix((i-1)*3+3)=1
osolve%w(i)=wstart+distzone*(wend-wstart)
!calculating temperature
select case(tempcase)
case (0)
osolve%kfixt(i)=1
osolve%temp(i)=1.d0
case (1)
osolve%kfixt(i)=1
osolve%temp(i)=osolve%x(i)*(1.0d0-temp2)+temp2
case (2)
osolve%kfixt(i)=1
osolve%temp(i)=distzone*(1.0d0-temp2)+temp2
case (3,4)
tempstart=(rx-dxy)/rx*(temp2-1.0d0)+1.0d0
tempend=1.0d0
osolve%kfixt(i)=1
osolve%temp(i)=tempstart+distzone*(tempend-tempstart)
case default
osolve%kfixt(i)=1
osolve%temp(i)=1.d0
end select
end if
else
!transition zone y-wise at both sides (symmetry)
if (osolve%x(i)>(rx+dxy)) then
!outside rotating cylinder x-wise
osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=vback !x-direction
osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0 !z-direction
select case(tempcase)
case (0,2,3,4)
osolve%kfixt(i)=1
osolve%temp(i)=1.d0
case (1)
osolve%kfixt(i)=1
osolve%temp(i)=osolve%x(i)*(1.0d0-temp2)+temp2
case default
osolve%kfixt(i)=1
osolve%temp(i)=1.d0
end select
elseif (osolve%x(i)<(rx-dxy)) then
!inside rotating cylinder x-wise
!values for start of transition zone
distmain=sqrt(osolve%x(i)**2+zmain**2)
theta=asin(osolve%x(i)/distmain)
ustart=distmain*omega*cos(theta)
uend=vback
wstart=-osolve%x(i)*omega
wend=0.d0
!set velocities
osolve%kfix((i-1)*3+1)=1
osolve%u(i)=ustart+distzone*(uend-ustart)
osolve%kfix((i-1)*3+3)=1
osolve%w(i)=wstart+distzone*(wend-wstart)
select case(tempcase)
case (0)
osolve%kfixt(i)=1
osolve%temp(i)=1.d0
case (1)
osolve%kfixt(i)=1
osolve%temp(i)=osolve%x(i)*(1.0d0-temp2)+temp2
case (2)
osolve%kfixt(i)=1
osolve%temp(i)=distzone*(1.0d0-temp2)+temp2
case (3,4)
tempstart=osolve%x(i)/rx*(temp2-1.0d0)+1.0d0
tempend=1.0d0
osolve%kfixt(i)=1
osolve%temp(i)=tempstart+distzone*(tempend-tempstart)
case default
osolve%kfixt(i)=1
osolve%temp(i)=1.d0
end select
else
!transition at corner
distmain=sqrt((rx-dxy)**2+zmain**2)
theta=asin((rx-dxy)/distmain)
!start velocities scale according to x-transition zone
ustart2=distmain*omega*cos(theta)
wstart2=-(rx-dxy)*omega
ustart=ustart2+(osolve%x(i)-(rx-dxy))/(2.d0*dxy)*(vback-ustart2)
wstart=wstart2+(osolve%x(i)-(rx-dxy))/(2.d0*dxy)*(-wstart2)
!y transition to (vback,0,0)
uend=vback
wend=0.d0
osolve%kfix((i-1)*3+1)=1
osolve%u(i)=ustart+distzone*(uend-ustart)
osolve%kfix((i-1)*3+3)=1
osolve%w(i)=wstart+distzone*(wend-wstart)
select case(tempcase)
case (0)
osolve%kfixt(i)=1
osolve%temp(i)=1.d0
case (1)
osolve%kfixt(i)=1
osolve%temp(i)=osolve%x(i)*(1.0d0-temp2)+temp2
case (2)
osolve%kfixt(i)=1
tempstart=temp2+(osolve%x(i)-(rx-dxy))/(2*dxy)*(1.0d0-temp2)
tempend=1.0d0
osolve%temp(i)=tempstart+distzone*(tempend-tempstart)
case (3,4)
tempstart2=(rx-dxy)/rx*(temp2-1.0d0)+1.0d0
tempstart=tempstart2+(osolve%x(i)-(rx-dxy))/(2*dxy)*(1.0d0-tempstart2)
tempend=1.0d0
osolve%kfixt(i)=1
osolve%temp(i)=tempstart+distzone*(tempend-tempstart)
case default
osolve%kfixt(i)=1
osolve%temp(i)=1.d0
end select
end if
end if
end if
!x=1 face with influx vback
if (osolve%x(i)>1.d0-eps) then
osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=vback
!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
!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
!case sphere
else
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)
!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
!strictflow: flow fixed to (vback, 0 0)
if (strictflow == 1) then
osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=vback
osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
end if
!rotating sphere cuts y=0 face
if (y0-0.5d0*wy < 0.d0) then
!calculate geometry for face
distedge = sqrt(max(0.0d0, rmain**2 - ((-y0)*rx/(0.5d0*wy))**2))
if (strictflow == 1) then
distzonerad = max(0.d0,min(1.d0,(distrot-distedge+dz)/(2*dz)))
else
if (distrot > distedge) then
distzonerad = 1.d0
else
distzonerad = 0.d0
end if
end if
!Setting velocities for inside and tranistion
!angle of rotation
theta=asin(osolve%x(i)/distrot)
!inside rotating cylinder
if (distzonerad < eps) then
!set x velocity
osolve%kfix((i-1)*3+1)=1
osolve%u(i)=distrot*omega*cos(theta)
!set z velocity
osolve%kfix((i-1)*3+3)=1
osolve%w(i)=-osolve%x(i)*omega
!in transition zone, only occurs for strictflow = 1
else if (distzonerad < 1.d0) then
xminedge = osolve%x(i)*(distedge-dz)/distrot
!set x velocity
ustart=(distedge-dz)*omega*cos(theta)
uend=vback
osolve%u(i)=ustart+distzonerad*(uend-ustart)
!set z velocity
wstart=-xminedge*omega
wend=0.d0
osolve%w(i)=wstart+distzonerad*(wend-wstart)
end if
end if
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
!strictflow: flow fixed to (vback, 0 0)
if (strictflow == 1) then
osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=vback
osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
end if
!rotating sphere cuts y=1 face
if (y0+0.5d0*wy > 1.d0) then
!calculate geometry for face
distedge = sqrt(max(0.0d0, rmain**2 - ((1.d0-y0)*rx/(0.5d0*wy))**2))
if (strictflow == 1) then
distzonerad = max(0.d0,min(1.d0,(distrot-distedge+dz)/(2*dz)))
else
if (distrot > distedge) then
distzonerad = 1.d0
else
distzonerad = 0.d0
end if
end if
!Setting velocities for inside and transition
!angle of rotation
theta=asin(osolve%x(i)/distrot)
!inside rotating cylinder
if (distzonerad < eps) then
!set x velocity
osolve%kfix((i-1)*3+1)=1
osolve%u(i)=distrot*omega*cos(theta)
!set z velocity
osolve%kfix((i-1)*3+3)=1
osolve%w(i)=-osolve%x(i)*omega
!in transition zone, only occurs for strictflow = 1
else if (distzonerad < 1.d0) then
xminedge = osolve%x(i)*(distedge-dz)/distrot
!set x velocity
ustart=(distedge-dz)*omega*cos(theta)
uend=vback
osolve%u(i)=ustart+distzonerad*(uend-ustart)
!set z velocity
wstart=-xminedge*omega
wend=0.d0
osolve%w(i)=wstart+distzonerad*(wend-wstart)
end if
end if
end if
!x=0 face
if (osolve%x(i)<eps) then
!set y and z velocity
osolve%kfix((i-1)*3+1)=1 !x-direction fixed
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
!outline of indenter
zedge = max(0.d0, zmain+sqrt(max(0.d0,rmain**2-((osolve%y(i)-y0)*rx/(0.5d0*wy))**2)))
zminedge = max(0.d0,zedge-dz)
zmaxedge = max(0.d0,zedge+dz)
!position relative to radial and y transition zones
! (0 inside, 1 outside, linear transition in zone)
distzonerad = max(0.d0,min(1.d0,(osolve%z(i)*params%vex-zminedge)/(zmaxedge-zminedge+eps)))
distzoney=min(1.d0,max(0.d0,(abs(osolve%y(i)-y0)-0.5d0*wy+dxy)/(2.d0*dxy)))
!outside rotating sphere radially
if (distzonerad > 1.d0-eps) then
ustart = vtop
uend = vback
osolve%u(i)=ustart+distzoney*(uend-ustart)
!within rotating sphere radially
elseif (distzonerad < eps) then
osolve%u(i)=omega*(osolve%z(i)*params%vex-zmain)
!inside transition zone radially
else
!uend depends on upper gradient (see distzonerad > 1)
ustart2=vtop
uend2=vback
uend=ustart2+distzoney*(uend2-ustart2)
!ustart is gradient from urot to uback
ustart2 = omega*(osolve%z(i)*params%vex-zmain)
uend2=vback
ustart=ustart2+distzoney*(uend2-ustart2)
osolve%u(i)=ustart+distzonerad*(uend-ustart)
end if
end if
!z=0 face
if (osolve%z(i)*params%vex < eps) then
!fix all velocities and temp, set v to 0
osolve%kfix((i-1)*3+1)=1
osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0 !y-direction
osolve%kfix((i-1)*3+3)=1
osolve%kfixt(i)=1
!calculating radial transition zone: radial belt of width 2dxy at indenter edge,
!circle origin is (x0, y0, 0)
!finding intersections with indenter edge radial direction
distz0plane = sqrt(osolve%x(i)**2+((osolve%y(i)-y0)*rx/(0.5d0*wy))**2)
if (distz0plane < eps) then
xedge = rx
yedge = y0
else
xedge = osolve%x(i)*rx/distz0plane
yedge = y0+(osolve%y(i)-y0)*rx/distz0plane
end if
!finding intersection with indenter edge for given x in y direction
yfix_xedge=sqrt(max(0.d0,rmain**2-((osolve%y(i)-y0)*rx/(0.5d0*wy))**2-zmain**2))
!start of transition zone radially
xminedge = xedge-dxy*xedge/sqrt(xedge**2+(yedge-y0)**2)
yminedge = yedge-dxy*(yedge-y0)/sqrt(xedge**2+(yedge-y0)**2)
!transition zone in radial >1 for outside, <0 for inside
distzonerad=(sqrt((osolve%y(i)-y0)**2+osolve%x(i)**2) - &
sqrt((yminedge-y0)**2+(xminedge)**2))/(2.d0*dxy)
!outside rotating sphere radially
if (distzonerad > 1.d0) then
!velocity
osolve%u(i) = vback
osolve%w(i) = 0.d0
!temperature
select case(tempcase)
case (1)
osolve%temp(i) = temp2+osolve%x(i)*(1.0d0-temp2)
case (2,3,4,5)
osolve%temp(i)=1.d0
case default
osolve%temp(i)=1.d0
end select
!within rotating sphere radially
elseif (distzonerad < 0.d0) then
distrot = sqrt(osolve%x(i)**2+zmain**2)
theta = asin(osolve%x(i)/distrot)
!set x velocity
osolve%u(i) = distrot*omega*cos(theta)
!set z velocity
osolve%w(i) = -osolve%x(i)*omega
select case(tempcase)
case (1)
osolve%temp(i) = temp2+osolve%x(i)*(1.0d0-temp2)
case (2)
osolve%temp(i) = temp2
case (3)
osolve%temp(i) = 1.d0+osolve%x(i)/rx*(temp2-1.0d0)
case (4)
osolve%temp(i) = 1.d0+osolve%x(i)/yfix_xedge*(temp2-1.0d0)
case (5)
osolve%temp(i) = 1.d0+distz0plane/rx*(temp2-1.0d0)
case default
osolve%temp(i)=1.d0
end select
!transition zone radially
else
!start value at inner edge of transition zone
distrot = sqrt(xminedge**2+zmain**2)
theta = asin(xminedge/distrot)
ustart = distrot*omega*cos(theta)
uend = vback
wstart = -xminedge*omega
wend = 0.d0
!setting velocities
osolve%u(i)=ustart+distzonerad*(uend-ustart)
osolve%w(i)=wstart+distzonerad*(wend-wstart)
!setting temperature
select case(tempcase)
case (1)
osolve%temp(i) = temp2+osolve%x(i)*(1.0d0-temp2)
case (2)
tempstart = temp2
tempend = 1.d0
osolve%temp(i)=tempstart+distzonerad*(tempend-tempstart)
case (3)
tempstart = 1.d0+xminedge/rx*(temp2-1.0d0)
tempend = 1.d0
osolve%temp(i)=tempstart+distzonerad*(tempend-tempstart)
case (4)
yfix_xedge = sqrt(max(0.d0,rmain**2-((yminedge-y0)*rx/(0.5d0*wy))**2-zmain**2))
tempstart = 1.d0+xminedge/yfix_xedge*(temp2-1.0d0)
tempend = 1.d0
osolve%temp(i)=tempstart+distzonerad*(tempend-tempstart)
case (5)
distz0plane = sqrt(xedge**2+((yedge-y0)*rx/(0.5d0*wy))**2)
tempstart = 1.d0+distz0plane/rx*(temp2-1.0d0)
tempend = 1.d0
osolve%temp(i)=tempstart+distzonerad*(tempend-tempstart)
case default
osolve%temp(i)=1.d0
end select
end if
end if
!x=1 face with influx vback
if (osolve%x(i)>1.d0-eps) then
osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=vback
if (strictflow == 1) then
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
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 if
end
!------------------------------------------------------------------------------|