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

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

            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)

              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,0,1) = c0

         diffus(i,0,2) = c0

         diffus(i,0,3) = c0

     I       ustar, bfsfc, hbl, stable, casea, diffus, kbl

      _RL diffus(imt,0:Nrp1,mdiff)

      _RL diffus(imt,0:Nrp1,mdiff)

         dvdzup = (diffus(i,kn-1,1) - diffus(i,kn  ,1)) / hwide(kn)

         dvdzdn = (diffus(i,kn  ,1) - diffus(i,kn+1,1)) / hwide(kn+1)

         dvdzup = (diffus(i,kn-1,2) - diffus(i,kn  ,2)) / hwide(kn)

         dvdzdn = (diffus(i,kn  ,2) - diffus(i,kn+1,2)) / hwide(kn+1)

         dvdzup = (diffus(i,kn-1,3) - diffus(i,kn  ,3)) / hwide(kn)

         dvdzdn = (diffus(i,kn  ,3) - diffus(i,kn+1,3)) / hwide(kn+1)

         visch  = diffus(i,kn,1) + viscp * delhat

         difsh  = diffus(i,kn,2) + difsp * delhat

         difth  = diffus(i,kn,3) + diftp * delhat

     I       dkm1, hbl, kbl, diffus, casea

      _RL diffus(imt,0:Nrp1,mdiff)

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

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

     O     , diffus, myThid )

             if(k.ge.kmtj(i))  diffus(i,k,md) = 0.0

     I       ustar, bfsfc, hbl, stable, casea, diffus, kbl

     I       dkm1, hbl, kbl, diffus, casea

     &              diffus(i,k,1) )

     &              diffus(i,k,2) )

     &              diffus(i,k,3) )

               diffus(i,k,1) = max ( blmc(i,k,1), viscArNr(1) )

               diffus(i,k,2) = max ( blmc(i,k,2), diffusKzS(i,Nr) )

               diffus(i,k,3) = max ( blmc(i,k,3), diffusKzT(i,Nr) )

     O     , diffus

     O       diffus,

      _RL diffus(imt,0:Nrp1,3)

      diffus(1,1,1) = 0.0

            diffus(i,ki,1) = 0.

            diffus(i,ki,2) = 0.

            diffus(i,ki,3) = 0.

               diffus(i,ki,1) = 0.

               diffus(i,ki,2) = 0.

               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))

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

     U     diffus(1,0,1),