integer i, kn, ki

      do ki = 1, Nr

         kkppkey = (ikppkey-1)*Nr + ki

            sig      = (-zgrid(ki) + 0.5 * hwide(ki)) / hbl(i)

            sig = (-zgrid(ki) + 0.5 * hwide(ki)) / hbl(i)

            blmc(i,ki,1) = hbl(i) * wm(i) * sig * (1. + sig * Gm)

            blmc(i,ki,2) = hbl(i) * ws(i) * sig * (1. + sig * Gs)

            blmc(i,ki,3) = hbl(i) * ws(i) * sig * (1. + sig * Gt)

            ghat(i,ki) = (1.-stable(i)) * cg / max(phepsi,tempVar)

      integer ki, i, md

         ki = kbl(i)-1

         if ((ki .ge. 1) .and. (ki .lt. Nr)) then

            delta = (hbl(i) + zgrid(ki)) / (zgrid(ki) - zgrid(ki+1))

               dkmp5         =      casea(i)  * diffus(i,ki,md) +

     1                         (1.- casea(i)) * blmc  (i,ki,md)

               blmc(i,ki,md) = (1.- delta)*diffus(i,ki,md)

            ghat(i,ki) = (1.- casea(i)) * ghat(i,ki)

      integer i, ki, kp1

      do ki = 1, Nr

            diffus(i,ki,1) = 0.

            diffus(i,ki,2) = 0.

            diffus(i,ki,3) = 0.

      do ki = 1, Nr

               diffus(i,ki,1) = 0.

               diffus(i,ki,2) = 0.

            elseif (ki      .GE. kmtj(i)) then

               diffus(i,ki,1) = diffus(i,ki-1,1)

               diffus(i,ki,2) = diffus(i,ki-1,2)

               diffus(i,ki,1) = dblocSm(i,ki) * (zgrid(ki)-zgrid(ki+1))

     &              / max( Shsq(i,ki), phepsi )

               diffus(i,ki,2) = dbloc(i,ki)   / (zgrid(ki)-zgrid(ki+1))

      do ki = 1, Nr

            Rig   = max ( diffus(i,ki,2) , BVSQcon )

            Rig  = max ( diffus(i,ki,1), c0 )

            kp1 = MIN(ki+1,Nr)

               diffus(i,ki,1) = viscArNr(1) + fcon*difmcon + fRi*difm0

               diffus(i,ki,2) = diffusKzS(i,kp1)+fcon*difscon+fRi*difs0

               diffus(i,ki,3) = diffusKzT(i,kp1)+fcon*diftcon+fRi*dift0

               diffus(i,ki,1) = viscArnr(1)

               diffus(i,ki,2) = diffusKzS(i,kp1)

               diffus(i,ki,3) = diffusKzT(i,kp1)

            diffus(i,ki,1) = viscArNr(1)

            diffus(i,ki,2) = diffusKzS(i,kp1)

            diffus(i,ki,3) = diffusKzT(i,kp1)