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

subroutine define_isostasy_bc(params,osolve,vo,zisodisp,firstcall,l,x0,spu,spv,&
           ldisp)

!------------------------------------------------------------------------------|
!((((((((((((((((  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

implicit none

type (parameters) params
type (octreesolve) osolve
type (void) vo
double precision zisodisp(2**params%levelmax_oct+1,2**params%levelmax_oct+1)
double precision x0(osolve%nnode),ldisp(osolve%nnode)
double precision l,spu,spv
logical firstcall

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

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

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

INCLUDE 'mpif.h'

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                                                                      ! Tiny 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+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)=(zisodisp(i+1,j)-zisodisp(i,j))/dxy                        ! Displacement between first and second nodes along x-axis
      elseif (i==nb+1) then
        zisoslx(i,j)=(zisodisp(i,j)-zisodisp(i-1,j))/dxy                        ! Displacement between second-to-last and last nodes along x-axis
      else
        zisoslx(i,j)=((zisodisp(i,j)-zisodisp(i-1,j))/dxy+&                     ! Average displacement of nodes on either side of current node, along x-axis
                     (zisodisp(i+1,j)-zisodisp(i,j))/dxy)/2.d0
      endif
      !if (j==1) then
      !  zisosly(i,j)=(zisodisp(i,j+1)-zisodisp(i,j))/dxy                        ! Displacement between first and second nodes along y-axis
      !elseif (j==nb+1) then
      !  zisosly(i,j)=(zisodisp(i,j)-zisodisp(i,j-1))/dxy                        ! Displacement between second-to-last and last nodes along y-axis
      !else
      !  zisosly(i,j)=((zisodisp(i,j)-zisodisp(i,j-1))/dxy+&                     ! Average displacement of nodes on either side of current node, along y-axis
      !               (zisodisp(i,j+1)-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)

    ! Modify velocities
    if (osolve%z(i).lt.eps) then
!      if (osolve%x(i).le.x0) then
!        if (spu.le.1.d0-eps) then                                                 ! No change in B/Cs
          if (osolve%u(i).gt.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
            osolve%u(i)=sqrt(vinit**2-osolve%w(i)**2.d0)                        ! New x-velocity
          endif
!        endif
!      else
!        if (spu.ge.eps) then
!          zdisp=osolve%u(i)*sin(atan(zsl))                                        ! Negative for negative slope
!          osolve%w(i)=zdisp-osolve%w(i)                                           ! New z-velocity
!          osolve%u(i)=sqrt(vin**2-osolve%w(i)**2.d0)                              ! New x-velocity
!        endif
!      endif
    endif
  enddo
  deallocate(zisoslx)
endif

! here we need to impose that the bc satisfy the bad faces constrains
!do ie=1,osolve%nface
!  do j=1,4
!    jp=1+mod(j,4)
!    i=osolve%iface(j+4,ie)
!    if (osolve%kfix((i-1)*3+1).eq.1) osolve%u(i)=(osolve%u(osolve%iface(j,ie))+&
!       osolve%u(osolve%iface(jp,ie)))/2.d0
!    if (osolve%kfix((i-1)*3+2).eq.1) osolve%v(i)=(osolve%v(osolve%iface(j,ie))+&
!       osolve%v(osolve%iface(jp,ie)))/2.d0
!    if (osolve%kfix((i-1)*3+3).eq.1) osolve%w(i)=(osolve%w(osolve%iface(j,ie))+&
!       osolve%w(osolve%iface(jp,ie)))/2.d0
!    if (osolve%kfixt(i).eq.1) osolve%temp(i)=(osolve%temp(osolve%iface(j,ie))+&
!       osolve%temp(osolve%iface(jp,ie)))/2.d0
!  enddo
!  i=osolve%iface(9,ie)
!  if (osolve%kfix((i-1)*3+1).eq.1) osolve%u(i)=&
!     sum(osolve%u(osolve%iface(1:4,ie)))/4.d0
!  if (osolve%kfix((i-1)*3+2).eq.1) osolve%v(i)=&
!     sum(osolve%v(osolve%iface(1:4,ie)))/4.d0
!  if (osolve%kfix((i-1)*3+3).eq.1) osolve%w(i)=&
!     sum(osolve%w(osolve%iface(1:4,ie)))/4.d0
!  if (osolve%kfixt(i).eq.1) osolve%temp(i)=&
!     sum(osolve%temp(osolve%iface(1:4,ie)))/4.d0
!enddo

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