!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              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
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
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

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
        
        !rotating sphere cuts y=0 face
		if (distmain <= rmain) then
			!angle of rotation
			theta=asin(osolve%x(i)/distrot)
			!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
		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
        
        !rotating cylinder cuts y=1 face
        if (distmain <= rmain) then
			distrot=sqrt(osolve%x(i)**2+(osolve%z(i)*params%vex-zmain)**2)
			!angle of rotation
			theta=asin(osolve%x(i)/distrot)
			!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
		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)
                    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 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 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
        !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

end if

end 
!------------------------------------------------------------------------------|
