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