Skip to content
Snippets Groups Projects
Commit 6e5ac32d authored by Dave Whipp's avatar Dave Whipp
Browse files

Attempted to generalize define_isostasy_bc.f90 to be able to be called by any... 

Attempted to generalize define_isostasy_bc.f90 to be able to be called by any defined velocity B/C option
parent 49b43a3e
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/define_isostasy_bc.f90
+ 64
− 86
View file @ 6e5ac32d
......@@ -16,8 +16,7 @@
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine define_isostasy_bc(params,osolve,zi,firstcall,l,x0,spu,spv,&
phi,ldisp)
subroutine define_isostasy_bc(params,osolve,bcdef)
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine )))))))))))))))))))))))))))))))))))))
......@@ -43,21 +42,15 @@ 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)
double precision l,spu,spv,phi
logical firstcall
type (bc_definition) bcdef
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
integer i,j,iproc,nproc,ierr,nb,xdisp,ydisp,xnow,ynow
!integer :: ie,jp
double precision eps,pi,zsl,dxy,zdisp,vinit,x1,x2
double precision,dimension(:,:),allocatable :: zisoslx
!double precision,dimension(:,:),allocatable :: zisosly
double precision,dimension(:,:),allocatable :: zisoslx,zisosly
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
......@@ -66,104 +59,89 @@ 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
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+zi%zisodisp(xdisp,ydisp)
enddo
else
allocate(zisoslx(nb+1,nb+1))!,zisosly(nb+1,nb+1))
!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)=(zi%zisodisp(i+1,j)-zi%zisodisp(i,j))/dxy ! Displacement between first and second nodes along x-axis
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)=(zi%zisodisp(i,j)-zi%zisodisp(i-1,j))/dxy ! Displacement between second-to-last and last nodes along x-axis
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
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
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
!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)
x1=x0(i)-l/tan(phi)
x2=x0(i)+l/tan(phi)
! 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
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
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)
deallocate(zisoslx,zisosly)
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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment