diff --git a/code/ysolve.f90 b/code/ysolve.f90
index 27c5b74..b17f903 100644
--- a/code/ysolve.f90
+++ b/code/ysolve.f90
@@ -86,6 +86,7 @@
       REAL :: c,yr,wrate,dely,sdr,sdd,sd,uu,onelw,dd
       INTEGER :: i,j,k,nd
       REAL, DIMENSION(1,nx) :: ux,dux,d2ux
+      REAL, DIMENSION(1,nx) :: duy
       REAL, DIMENSION(10,nx) :: sav
         
       sav=0.
@@ -103,12 +104,14 @@
        IF (ux(1,i).gt.1.) ux(1,i)=1.                                  
       ENDDO
           
+      dm = diff * ((1.0 + bc * (1.0 - y1(:,1,1))) ** 0.76)
+
       nd=1
       CALL dfnonp(nx,hx,ux(1,:),dux(1,:),nd,1)  
+      CALL dfnonp(nx,hx,dm,duy(1,:),nd,1)
       nd=1 
       CALL d2fnonp(nx,hx,ux(1,:),d2ux(1,:),nd,1) 
 
-      dm = diff * (1.0 + (rvis - 1.0) * (1.0 - y1(:,1,1)))
 
       DO i=1,nx
        yr=ux(1,i)
@@ -131,14 +134,14 @@
 
        dely=ABS(dux(1,i))
        sdr=wrate/dely
-       sdd=-dm(i)*d2ux(1,i)/dely
+       sdd=-(dm(i)*d2ux(1,i)+duy(1,i)*dux(1,i))/dely
        IF (dely.eq.0.) THEN
         sdr=0.; sdd=0.
        ENDIF
        sd=sdr+sdd
        uu=u(1,1,1)
        onelw=(-dux(1,i))/c
-       dd=-dm(i)*d2ux(1,i)
+       dd=-(dm(i)*d2ux(1,i)+duy(1,i)*dux(1,i))
        if (c.eq.0.) onelw=0.
 
        sav(1,i)=sav(1,i)+yr
@@ -373,6 +376,7 @@
       REAL, INTENT(IN),DIMENSION(:,:,:) :: r1
       REAL, INTENT(OUT),DIMENSION(:,:,:) :: f
       REAL, DIMENSION(3,nx) :: ux,dux,d2ux
+      REAL, DIMENSION(1,nx) :: duy
       INTEGER :: i,j,k
       REAL :: wrate,Ly,Dy
 
@@ -385,10 +389,11 @@
          ux(3,i)=r1(i,j,k)                              ! Y
         ENDDO
 
-        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 * ((1.0 + bc * (1.0 - r1(:,1,1))) ** 0.76)
 
-      dm = diff * (1.0 + (rvis - 1.0) * (1.0 - r1(:,1,1)))
+        CALL dfnonp(nx,hx,ux(1:3,:),dux(1:3,:),3,1)
+        CALL dfnonp(nx,hx,dm,duy(1,:),1,1)
+        CALL d2fnonp(nx,hx,ux(3:3,:),d2ux(1,:),1,1)
 
         DO i=1,nx
          wrate=pre*ux(3,i)*exp(-ac/(1.+bc*(1.-ux(3,i))))   !wrate
@@ -406,7 +411,7 @@
 ! -0.5*( d(u*Y)/dx + u*dY/dx + Y*du/dx ) + D*d2Y/d2x
          f(i,j,k)=-0.5*( dux(1,i) + ux(2,i)*dux(3,i) + &
           ux(3,i)*dux(2,i) ) &
-                  + dm(i)*d2ux(1,i) - wrate
+                  + dm(i)*d2ux(1,i) + duy(1,i)*dux(3,i) - wrate
          IF (i.eq.nx) THEN
           Ly=ux(2,nx)*dux(3,nx)                           ! Ly = u*dYr/dx
           f(nx,1,1)=-wrate