Modified to only apply B/Cs on either side of s-point (x0)

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)