      SUBROUTINE fkdat

      IMPLICIT NONE

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

      INTEGER iseed
      REAL*4 ran
      REAL ranf

      REAL alpha,beta,preexp,incrho,rhoi

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

      iseed=1.0

      DO i=1,imax
         DO j=1,jmax
            DO k=1,kmax

               unew(i,j,k) = 0.0
               vnew(i,j,k) = 0.0
               wnew(i,j,k) = 0.1

c..............SGI machine.........

c               vnew(i,j,k) = -0.05 + .1*ran(iseed)
c               unew(i,j,k) = -0.05 + .1*ran(iseed)

c..............Cray...............

c               vnew(i,j,k) = -0.05 + .1*ranf()
c               unew(i,j,k) = -0.05 + .1*ranf()
c               wnew(i,j,k) = -0.05 + .1*ranf()

               tnew(i,j,k) = 0.0
               told(i,j,k) = 0.0
               ynew(i,j,k) = 0.0

               munew(i,j,k)  = 1.0
               lmbnew(i,j,k) = 1.0
               difnew(i,j,k) = 1.0

               mut(i,j,k)    = 0.0

               dudtn(i,j,k)  = 0.0
               dudtnm(i,j,k) = 0.0
               dvdtn(i,j,k)  = 0.0
               dvdtnm(i,j,k) = 0.0
               dwdtn(i,j,k)  = 0.0
               dwdtnm(i,j,k) = 0.0
               dtdtn(i,j,k)  = 0.0
               dtdtnm(i,j,k) = 0.0
               dydtn(i,j,k)  = 0.0
               dydtnm(i,j,k) = 0.0

            END DO
         END DO
      END DO

      DO k=0,k1
         DO j=0,j1
            DO i=0,i1
               rhonew(i,j,k) = 1.0/((incrho*alpha/(1.0-incrho*alpha))
     $              *tnew(i,j,k)+rhoi)
               rhold1(i,j,k) = rhonew(i,j,k)
               rhold2(i,j,k) = rhonew(i,j,k)
            END DO
         END DO
      END DO

      DO k=1,kmax
         DO j=1,jmax
            DO i=1,imax

               umnew(i,j,k) = rhonew(i,j,k)*unew(i,j,k)
               vmnew(i,j,k) = rhonew(i,j,k)*vnew(i,j,k)
               wmnew(i,j,k) = rhonew(i,j,k)*wnew(i,j,k)
               ymnew(i,j,k) = rhonew(i,j,k)*ynew(i,j,k)

               umold(i,j,k) = rhonew(i,j,k)*unew(i,j,k)
               vmold(i,j,k) = rhonew(i,j,k)*vnew(i,j,k)
               wmold(i,j,k) = rhonew(i,j,k)*wnew(i,j,k)
               ymold(i,j,k) = rhonew(i,j,k)*ynew(i,j,k)
               
            END DO
         END DO
      END DO

      RETURN
      END
