Newer
Older
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! DEFINE_ISOSTASY_BC OCT. 2009 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine define_isostasy_bc(params,osolve,zi,firstcall,l,x0,spu,spv,&
Dave Whipp
committed
phi,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
!use mpi
include 'mpif.h'
type (parameters) params
type (octreesolve) osolve
type (ziso) zi
!double precision zisodisp(2**params%levelmax_oct+1,2**params%levelmax_oct+1)
double precision x0(osolve%nnode),ldisp(osolve%nnode)
Dave Whipp
committed
double precision l,spu,spv,phi
logical firstcall
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
integer i,j,iproc,nproc,ierr,nb,xdisp,ydisp,xnow,ynow
!integer :: ie,jp
double precision,dimension(:,:),allocatable :: zisoslx
!double precision,dimension(:,:),allocatable :: 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 ! 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+zi%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)=(zi%zisodisp(i+1,j)-zi%zisodisp(i,j))/dxy ! Displacement between first and second nodes along x-axis
zisoslx(i,j)=(zi%zisodisp(i,j)-zi%zisodisp(i-1,j))/dxy ! Displacement between second-to-last and last nodes along x-axis
zisoslx(i,j)=((zi%zisodisp(i,j)-zi%zisodisp(i-1,j))/dxy+& ! Average displacement of nodes on either side of current node, along x-axis
(zi%zisodisp(i+1,j)-zi%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)
Dave Whipp
committed
x1=x0(i)-l/tan(phi)
x2=x0(i)+l/tan(phi)
! Modify velocities
if (osolve%z(i).lt.eps) then
Dave Whipp
committed
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
Dave Whipp
committed
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
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|