diff --git a/code/ysolve.f90 b/code/ysolve.f90
index b17f903..c73bc27 100644
--- a/code/ysolve.f90
+++ b/code/ysolve.f90
@@ -104,7 +104,7 @@
        IF (ux(1,i).gt.1.) ux(1,i)=1.                                  
       ENDDO
           
-      dm = diff * ((1.0 + bc * (1.0 - y1(:,1,1))) ** 0.76)
+      dm = (/ (diff * diffusivity_T076(1.-y1(i,1,1)), i=1,nx) /)
 
       nd=1
       CALL dfnonp(nx,hx,ux(1,:),dux(1,:),nd,1)  
@@ -389,7 +389,7 @@
          ux(3,i)=r1(i,j,k)                              ! Y
         ENDDO
 
-      dm = diff * ((1.0 + bc * (1.0 - r1(:,1,1))) ** 0.76)
+      dm = (/ (diff * diffusivity_T076(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)
@@ -470,5 +470,21 @@
       f(nx,1,1)=f(nx,1,1)-Dy
 
      END SUBROUTINE fns
+!------------------------------------------------------------------------
+     REAL FUNCTION sigmoid(x)
+         REAL, INTENT(IN) :: x
+         sigmoid = 1. / (1 + exp(-x))
+     END FUNCTION sigmoid
+!------------------------------------------------------------------------
+     REAL FUNCTION diffusivity_sigmoid(c)
+         REAL, INTENT(IN) :: c
+         REAL :: a = 200.
+         diffusivity_sigmoid = 1. + (rvis - 1.) * sigmoid(a*(c-cstar))
+     END FUNCTION diffusivity_sigmoid
+!------------------------------------------------------------------------
+     REAL FUNCTION diffusivity_T076(c)
+         REAL, INTENT(IN) :: c
+         diffusivity_T076 = ((1.0 + bc * c) ** 0.76)
+     END FUNCTION diffusivity_T076
 !------------------------------------------------------------------------
     END MODULE ysolve