diff --git a/code/ysolve.f90 b/code/ysolve.f90 index 3efdf14..ae3f5bf 100644 --- a/code/ysolve.f90 +++ b/code/ysolve.f90 @@ -108,6 +108,9 @@ nd=1 CALL d2fnonp(nx,hx,ux(1,:),d2ux(1,:),nd,1) + dm = diff + where (y1(:,1,1)<(1-cstar)) dm = diff * rvis + DO i=1,nx yr=ux(1,i) c=1.-yr @@ -386,6 +389,9 @@ CALL dfnonp(nx,hx,ux(1:3,:),dux(1:3,:),3,1) CALL d2fnonp(nx,hx,ux(3:3,:),d2ux(1,:),1,1) + dm = diff + where (r1(:,1,1)<(1.0-cstar)) dm = diff * rvis + DO i=1,nx wrate=pre*ux(3,i)*exp(-ac/(1.+bc*(1.-ux(3,i)))) !wrate ! IF ((1.-ux(3,i)).le.c_cut) THEN