!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!            DEFINE_ISOSTASY_BC   OCT. 2009                                    |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

subroutine define_isostasy_bc(params,osolve,bcdef)

!------------------------------------------------------------------------------|
!((((((((((((((((  Purpose of the routine  )))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! New subroutine to apply the modified isostasy velocity boundary conditions
!  - This routine uses the basal displacement array defined in isostasy.f90 to
!    modify the basal velocity boundary conditions to have the incoming flux
!    velocity tangent to the moho. This might not be exactly correct, but works
!    for now.
!
! dwhipp (11/09)
!
!------------------------------------------------------------------------------|
!((((((((((((((((  declaration of the subroutine arguments  ))))))))))))))))))))
!------------------------------------------------------------------------------|

use definitions
!use mpi

implicit none

include 'mpif.h'

type (parameters) params
type (octreesolve) osolve
type (bc_definition) bcdef

!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|

integer i,j,iproc,nproc,ierr,nb,xdisp,ydisp,xnow,ynow
double precision eps,pi,zsl,dxy,zdisp,vinit,x1,x2
double precision,dimension(:,:),allocatable :: zisoslx,zisosly

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)

! General variable definition
eps=1.d-10                                                                      ! Very small number
pi=atan(1.d0)*4.d0                                                              ! Pi
nb=2**params%levelmax_oct                                                       ! Number of nodes along one side of model
dxy=1./real(nb)                                                                 ! Non-dimensional node spacing

!if (firstcall) then
!  do i=1,osolve%nnode
!    ! Shift s-line depth
!    xdisp=nint(x0(i)*nb)+1
!    ydisp=nint(osolve%y(i)*nb)+1
!    ldisp(i)=l+bcdef%zisodisp(xdisp,ydisp)
!  enddo
!else
  allocate(zisoslx(nb+1,nb+1)),zisosly(nb+1,nb+1))  
  do j=1,nb+1
    do i=1,nb+1
      if (i==1) then
        zisoslx(i,j)=(bcdef%zisodisp(i+1,j)-bcdef%zisodisp(i,j))/dxy            ! Displacement between first and second nodes along x-axis
      elseif (i==nb+1) then
        zisoslx(i,j)=(bcdef%zisodisp(i,j)-bcdef%zisodisp(i-1,j))/dxy            ! Displacement between second-to-last and last nodes along x-axis
      else
        zisoslx(i,j)=((bcdef%zisodisp(i,j)-bcdef%zisodisp(i-1,j))/dxy+&         ! Average displacement of nodes on either side of current node, along x-axis
                     (bcdef%zisodisp(i+1,j)-bcdef%zisodisp(i,j))/dxy)/2.d0
      endif
      if (j==1) then
        zisosly(i,j)=(bcdef%zisodisp(i,j+1)-bcdef%zisodisp(i,j))/dxy            ! Displacement between first and second nodes along y-axis
      elseif (j==nb+1) then
        zisosly(i,j)=(bcdef%zisodisp(i,j)-bcdef%zisodisp(i,j-1))/dxy            ! Displacement between second-to-last and last nodes along y-axis
      else
        zisosly(i,j)=((bcdef%zisodisp(i,j)-bcdef%zisodisp(i,j-1))/dxy+&         ! Average displacement of nodes on either side of current node, along y-axis
                     (bcdef%zisodisp(i,j+1)-bcdef%zisodisp(i,j))/dxy)/2.d0
      endif
    enddo
  enddo

  do i=1,osolve%nnode
    xnow=nint(osolve%x(i)*nb)+1
    ynow=nint(osolve%y(i)*nb)+1
!    zsl=zisoslx(xnow,ynow)
!    x1=x0(i)-l/tan(phi)
!    x2=x0(i)+l/tan(phi)

    ! Modify velocities
    if (osolve%z(i).lt.eps) then
      uv=sqrt(osolve%u(i)**2.d0+osolve%v(i)**2.d0)
      uvw=sqrt(osolve%u(i)**2.d0+osolve%v(i)**2.d0+osolve%w(i)**2.d0)
      thetax=atan(zisoslx(xnow,ynow))
      thetay=atan(zisosly(xnow,ynow))
      newu=
      osolve%u(i)=osolve%u(i)+(1-cos(thetax))
      osolve%v(i)=osolve%v(i)+(1-cos(thetay))
      ! PLUS OR MINUS HERE???
      osolve%w(i)=osolve%w(i)+((osolve%u(i)/uv)*sin(thetax)+(osolve%v(i)/uv)*  &
                  sin(thetay))

!      if (osolve%x(i).le.x1) then
!        if (spu.le.1.d0-eps) then                                               ! No change in B/Cs
!          zdisp=osolve%u(i)*sin(atan(zsl))                                      ! Negative for negative slope
!          vinit=(osolve%u(i)**2.d0+osolve%w(i)**2.d0)**0.5d0                    ! Incoming velocity vector magnitude
!          osolve%w(i)=osolve%w(i)+zdisp                                         ! New z-velocity
!          if (osolve%u(i).ge.0.d0) then
!            osolve%u(i)=sqrt(vinit**2-osolve%w(i)**2.d0)                        ! New x-velocity
!          else
!            osolve%u(i)=-sqrt(vinit**2-osolve%w(i)**2.d0)                       ! New x-velocity
!          endif
!        endif
!      elseif (osolve%x(i).ge.x2) then
!        if (spu.ge.eps) then
!          zdisp=osolve%u(i)*sin(atan(zsl))                                      ! Negative for negative slope
!          vinit=(osolve%u(i)**2.d0+osolve%w(i)**2.d0)**0.5d0                    ! Incoming velocity vector magnitude
!          osolve%w(i)=osolve%w(i)+zdisp                                         ! New z-velocity
!          if (osolve%u(i).ge.0.d0) then
!            osolve%u(i)=sqrt(vinit**2-osolve%w(i)**2.d0)                        ! New x-velocity
!          else
!            osolve%u(i)=-sqrt(vinit**2-osolve%w(i)**2.d0)                       ! New x-velocity
!          endif
!        endif
!      endif
    endif
  enddo
  deallocate(zisoslx,zisosly)
endif

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