From ea6f6644d1c698e45864903deb4c6e0c93695e09 Mon Sep 17 00:00:00 2001
From: Dave Whipp <dwhipp@dal.ca>
Date: Tue, 16 Feb 2010 18:38:24 +0000
Subject: [PATCH] Modified isostasy B/Cs to only flux extra material out before
the mantle subducts
---
src/define_isostasy_bc.f90 | 10 ++++++----
1 file changed, 6 insertions(+), 4 deletions(-)
diff --git a/src/define_isostasy_bc.f90 b/src/define_isostasy_bc.f90
index d964ec25..df17ee18 100644
--- a/src/define_isostasy_bc.f90
+++ b/src/define_isostasy_bc.f90
@@ -17,7 +17,7 @@
!------------------------------------------------------------------------------|
subroutine define_isostasy_bc(params,osolve,vo,zi,firstcall,l,x0,spu,spv,&
- ldisp)
+ phi,ldisp)
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine )))))))))))))))))))))))))))))))))))))
@@ -44,7 +44,7 @@ type (void) vo
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
+double precision l,spu,spv,phi
logical firstcall
!------------------------------------------------------------------------------|
@@ -103,10 +103,12 @@ else
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
- if (osolve%x(i).le.x0(i)) 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
@@ -117,7 +119,7 @@ else
osolve%u(i)=-sqrt(vinit**2-osolve%w(i)**2.d0) ! New x-velocity
endif
endif
- else
+ 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
--
GitLab