diff --git a/code/ysolve.f90 b/code/ysolve.f90
index 8d835ff..cbcff2b 100644
--- a/code/ysolve.f90
+++ b/code/ysolve.f90
@@ -134,7 +134,7 @@
        IF (ux(1,i).gt.1.) ux(1,i)=1.                                  
       ENDDO
           
-      dm = (/ (diff * diffusivity_T076(1.-y1(i,1,1)), i=1,nx) /)
+      dm = (/ (diff * diffusivity_sutherland(1.-y1(i,1,1)), i=1,nx) /)
 
       nd=1
       CALL dfnonp(nx,hx,ux(1,:),dux(1,:),nd,1)  
@@ -335,14 +335,6 @@
         WRITE(otape,*) cdum,prof_wr
         READ(itape,*) cdum,rvis
         WRITE(otape,*) cdum,rvis
-        READ(itape,*) cdum,cstar
-        WRITE(otape,*) cdum,cstar
-        READ(itape,*) cdum,sigw
-        WRITE(otape,*) cdum,sigw
-
-       DO i = 0,100
-        WRITE(*,*) diffusivity_T076(i / 100.)
-       END DO
 
        IF (.not.read_stdin) THEN
          CLOSE(itape)
@@ -432,7 +424,7 @@
          ux(3,i)=r1(i,j,k)                              ! Y
         ENDDO
 
-      dm = (/ (diff * diffusivity_T076(1.-r1(i,1,1)), i=1,nx) /)
+      dm = (/ (diff * diffusivity_sutherland(1.-r1(i,1,1)), i=1,nx) /)
 
         CALL dfnonp(nx,hx,ux(1:3,:),dux(1:3,:),3,1)
         CALL dfnonp(nx,hx,dm,duy(1,:),1,1)
@@ -485,5 +477,19 @@
          REAL, INTENT(IN) :: c
          diffusivity_T076 = ((1.0 + bc * c) ** 0.76)
      END FUNCTION diffusivity_T076
+!------------------------------------------------------------------------
+     REAL FUNCTION diffusivity_sutherland(c)
+         REAL, INTENT(IN) :: c
+         REAL :: theta, As, Ts, T0, T1
+         T0 = 1.0
+         T1 = (1.0 + bc)
+
+         Ts = (rvis*T1 - (T1**(3./2.))) / (T1**(3./2.) - rvis)
+         As = (T0 + Ts)
+
+         theta = (1.0 + bc * c)
+
+         diffusivity_sutherland = As * sqrt(theta) / (1. + Ts/theta)
+     END FUNCTION diffusivity_sutherland
 !------------------------------------------------------------------------
     END MODULE ysolve