diff --git a/src/define_isostasy_bc.f90 b/src/define_isostasy_bc.f90 index d964ec25eb09f46e951353f0d044c3bbd8e950f7..df17ee1876f97d2782d10224e5d5c8fbd7fc10fe 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