Modified isostasy B/Cs to only flux extra material out before the mantle subducts

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