      SUBROUTINE diffu2(ib,ie,jb,je,kb,ke)

      IMPLICIT NONE

      INCLUDE 'param.txt'
      INCLUDE 'common.txt'

      INTEGER  ib,ie,jb,je,kb,ke

c********************************************************************
c
c     Diffus calculates the diffusion for all variables over the
c     indicated domain. The indices ib,ie,jb,je,kb and ke ("i,j,k" for
c     the direction and "b" and "e" for begin and end) determine the
c     points for which the diffusion is calculated.
c
c     In this numerical model, we substituted the Euler-forward scheme
c     by a partly explicit and partly implicit method. In a fully
c     explicit scheme (like Euler-scheme) stability conditions limit
c     the maximum admissible timestep in a very severe way. Therefore,
c     we treat the diffusion in r- and z-direction explicitly in time,
c     whereas we treat the diffusion in phi-direction implicitly in time
c     to avoid stability problems.
c
c*********************************************************************


      CALL diffua(dudtnm,unew,vnew,wnew,munew,ru,rp,dphi,zw,zp,i1,j1,k1
     $     ,ib,ie,jb,je,kb,ke)

      CALL diffva(dvdtnm,unew,vnew,wnew,munew,ru,rp,dr,dphi,zw,zp,i1,j1
     $     ,k1,ib,ie,jb,je,kb,ke)

      CALL diffwa(dwdtnm,unew,vnew,wnew,munew,ru,rp,dphi,zw,zp,i1,j1,k1
     $     ,ib,ie,jb,je,kb,ke)

      CALL difftb(dtdtnm,unew,vnew,wnew,lmbnew,tnew,ynew,rhonew,ru,rp
     $     ,dphi,zw,zp,i1,j1,k1,ib,ie,jb,je,kb,ke,time)

      CALL diffyb(dydtnm,unew,vnew,wnew,difnew,tnew,ynew,rhonew,ru,rp
     $     ,dphi,zw,zp,i1,j1,k1,ib,ie,jb,je,kb,ke,time)

      RETURN
      END 

c--------------------------------------------------------------------------
c--------------------------------------------------------------------------
      SUBROUTINE difftb(dtdt,u,v,w,lmb,t,y,rho,ru,rp,dphi,zw,zp,i1,j1,k1
     $     ,ib,ie,jb,je,kb,ke,time)

      IMPLICIT NONE

      INCLUDE 'param2.txt'

      INTEGER i1,j1,k1,ib,ie,jb,je,kb,ke
      INTEGER i,j,k,im,ip,jm,jp,km,kp

      REAL t(0:i1,0:j1,0:k1),rho(0:i1,0:j1,0:k1)
      REAL dtdt(0:i1,0:j1,0:k1),u(0:i1,0:j1,0:k1)
      REAL v(0:i1,0:j1,0:k1),w(0:i1,0:j1,0:k1)
      REAL lmb(0:i1,0:j1,0:k1),y(0:i1,0:j1,0:k1)

      REAL ru(0:i1),rp(0:i1),dphi,zw(0:k1),zp(0:k1)

      REAL Rei,Pri

      REAL time

      REAL Re,Pr,Le

      REAL alpha,beta,preexp,incrho,rhoi

      COMMON /nonnum/ Re,Pr,Le

      COMMON /heatrl/ alpha,beta,preexp,incrho,rhoi
c----------------------------------------------------------------------

      DO k=kb,ke

         kp=k+1
         km=k-1

         DO j=jb,je

            jp=j+1
            jm=j-1

            DO i=ib,ie

               ip=i+1
               im=i-1

               dtdt(i,j,k)=dtdt(i,j,k) + (1.0/rho(i,j,k))*(
     1              (ru(i)*(lmb(i,j,k)+lmb(ip,j,k))*(t(ip,j,k)-t(i,j,k
     1              ))/(rp(ip)-rp(i)) - ru(im)*(lmb(i,j,k)+lmb(im,j,k)
     1              )*(t(i,j,k)-t(im,j,k))/(rp(i)-rp(im)))/(2.0*rp(i)
     1              *(ru(i)-ru(im)))
     2              +((lmb(i,jp,k)+lmb(i,j,k))*(t(i,jp,k)-t(i,j,k)) -
     2              (lmb(i,j,k)+lmb(i,jm,k))*(t(i,j,k)-t(i,jm,k)))/(2
     2              .0*rp(i)**2*dphi**2)
     3              +((lmb(i,j,kp)+lmb(i,j,k))*(t(i,j,kp)-t(i,j,k))
     3              /(zp(kp)-zp(k)) - (lmb(i,j,k)+lmb(i,j,km))*(t(i,j
     3              ,k)-t(i,j,km))/(zp(k)-zp(km)))/(2.0*(zw(k)-zw(km))))
     4                    +
     4              preexp*y(i,j,k)*EXP((-beta*(1.0-t(i,j
     4              ,k)))/(1-alpha*(1.0-t(i,j,k))))


            END DO

         END DO

      END DO

      RETURN
      END 

c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
      SUBROUTINE diffyb(dydt,u,v,w,dif,t,y,rho,ru,rp,dphi,zw,zp,i1,j1,k1
     $     ,ib,ie,jb,je,kb,ke,time)

      IMPLICIT NONE

      INCLUDE 'param2.txt'

      INTEGER i1,j1,k1,ib,ie,jb,je,kb,ke
      INTEGER i,j,k,im,ip,jm,jp,km,kp

      REAL y(0:i1,0:j1,0:k1),rho(0:i1,0:j1,0:k1)
      REAL dydt(0:i1,0:j1,0:k1),u(0:i1,0:j1,0:k1)
      REAL v(0:i1,0:j1,0:k1),w(0:i1,0:j1,0:k1)
      REAL dif(0:i1,0:j1,0:k1),t(0:i1,0:j1,0:k1)

      REAL ru(0:i1),rp(0:i1),dphi,zw(0:k1),zp(0:k1)

      REAL Rei,Pri,Lei

      REAL time

      REAL Re,Pr,Le

      REAL alpha,beta,preexp,incrho,rhoi

      COMMON /nonnum/ Re,Pr,Le

      COMMON /heatrl/ alpha,beta,preexp,incrho,rhoi
c----------------------------------------------------------------------

      DO k=kb,ke

         kp=k+1
         km=k-1

         DO j=jb,je

            jp=j+1
            jm=j-1

            DO i=ib,ie

               ip=i+1
               im=i-1

               dydt(i,j,k)=dydt(i,j,k) + (
     1              (ru(i)*(dif(i,j,k)+dif(ip,j,k))*(y(ip,j,k)-y(i,j,k
     1              ))/(rp(ip)-rp(i)) - ru(im)*(dif(i,j,k)+dif(im,j,k)
     1              )*(y(i,j,k)-y(im,j,k))/(rp(i)-rp(im)))/(2.0*rp(i)
     1              *(ru(i)-ru(im)))
     2              +((dif(i,jp,k)+dif(i,j,k))*(y(i,jp,k)-y(i,j,k)) -
     2              (dif(i,j,k)+dif(i,jm,k))*(y(i,j,k)-y(i,jm,k)))/(2
     2              .0*rp(i)**2*dphi**2)
     3              +((dif(i,j,kp)+dif(i,j,k))*(y(i,j,kp)-y(i,j,k))
     3              /(zp(kp)-zp(k)) - (dif(i,j,k)+dif(i,j,km))*(y(i,j
     3              ,k)-y(i,j,km))/(zp(k)-zp(km)))/(2.0*(zw(k)-zw(km))))
     4                    -
     4               preexp*rho(i,j,k)*y(i,j,k)*EXP((-beta*(1.0-t(i,j
     4              ,k)))/(1-alpha*(1.0-t(i,j,k))))

            END DO

         END DO

      END DO

      RETURN
      END 

