!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! 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(i)) 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
else
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)
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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|