From 2c9074a9c4ec34ce5ccce2b8ef2bdf722f224dd3 Mon Sep 17 00:00:00 2001
From: Dave Whipp <dwhipp@dal.ca>
Date: Wed, 23 Dec 2009 18:34:09 +0000
Subject: [PATCH] Modified to only apply B/Cs on either side of s-point (x0)

---
 src/define_isostasy_bc.f90 | 35 +++++++++++++++++++++--------------
 1 file changed, 21 insertions(+), 14 deletions(-)

diff --git a/src/define_isostasy_bc.f90 b/src/define_isostasy_bc.f90
index 65b26e4f..3f8ba11f 100644
--- a/src/define_isostasy_bc.f90
+++ b/src/define_isostasy_bc.f90
@@ -105,22 +105,29 @@ else
 
     ! Modify velocities
     if (osolve%z(i).lt.eps) then
-!      if (osolve%x(i).le.x0) then
-!        if (spu.le.1.d0-eps) then                                                 ! No change in B/Cs
-          if (osolve%u(i).gt.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%x(i).le.x0) 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
-!          osolve%w(i)=zdisp-osolve%w(i)                                           ! New z-velocity
-!          osolve%u(i)=sqrt(vin**2-osolve%w(i)**2.d0)                              ! New x-velocity
-!        endif
-!      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)
-- 
GitLab