C***********************************************************
C***********************************************************
C
      SUBROUTINE  NLDIFA
     I   (   QADV, idvar,
     I       dt, dx,
     I       IDF,
     O       DTQDNL,
     W       QLAP, CDQNLX, CDQNLZ,
     B       zgrid, DENS0, FDZ,
     D       NX, NZ, NXB, NZB, NROT )
C
C***********************************************************
      IMPLICIT REAL*8 ( A-H, O-Z )
C***********************************************************
      parameter ( nzmax = 1002, nvar = 10 )
c
      dimension sclap( -2:nzmax, nvar )
      save sclap
C
      DIMENSION
     \     QADV  ( -2:NXB, -2:NZB, NROT ),
     \     DTQDNL( -2:NXB, -2:NZB ),
     \     QLAP  ( -2:NXB, -2:NZB ),
     \     CDQNLX( -2:NXB, -2:NZB ),
     \     CDQNLZ( -2:NXB, -2:NZB ),
     \     DENS0 (-2:NZB), FDZ(-2:NZB),
     \     ZGRID (-2:NZB)
C
C=====================================================================
C     COMPUTATIONAL DIFFUSION TERM ( ARTIFICIAL SMOOTHING )
C --------------------------------------------------------------------
C
C *** Calculation of pseudo Laplacian ***
C
      DO 11 IZ = 0, NZ-1
      DO 11 IX = 1, NX
         QLAP(IX,IZ) =
     \ (3.*ABS(QADV(IX+1,IZ,IDF)+QADV(IX-1,IZ,IDF)-2.*QADV(IX,IZ,IDF))
     \ +ABS(QADV(IX,IZ+1,IDF)+QADV(IX,IZ-1,IDF)-2.*QADV(IX,IZ,IDF)) )
     \ / sclap(iz,idvar)

   11 CONTINUE
C
      CALL BOUND ( QLAP, NX, NZ, NXB, NZB, 1, 1, 1, 1 )
C
C *** Calculation of Diffusivity ***
C
      CD00 =    ( dx**2 / dt ) * 0.01D0 
      CDMAX =   ( dx**2 / dt ) / 5.0d0
C
      DO 21 IZ = 0, NZ
      DO 21 IX = 1, NX
         CDQNLX(IX,IZ) = MIN ( CDMAX, (
     \      CD00 * ( ( QLAP(IX,IZ) + QLAP(IX-1,IZ) ) / 2. )  ) )
         CDQNLZ(IX,IZ) = MIN ( CDMAX, (
     \      CD00 * ( ( QLAP(IX,IZ) + QLAP(IX,IZ-1) ) / 2. )  ) )
   21 CONTINUE
C
      CALL BOUND( CDQNLX, NX, NZ, NXB, NZB, 1, 1, 1, 1 )
      CALL BOUND( CDQNLZ, NX, NZ, NXB, NZB, 1, 1, 1, 1 )
C
C *** Time derivatives due to Non-linear Diffusion ***
C
      DO 31 IZ = 0, NZ-1
      DO 31 IX = 1, NX

         DQDINX = 
     \   - CDQNLX(IX  ,IZ) * ( QADV(IX  ,IZ,IDF) - QADV(IX-1,IZ,IDF) )
     \   + CDQNLX(IX+1,IZ) * ( QADV(IX+1,IZ,IDF) - QADV(IX  ,IZ,IDF) )

C
         DQDINZ =
     \    (  ( DENS0(IZ+1) + DENS0(IZ) ) / 2.0D0 
     \         * CDQNLZ(IX,IZ+1) 
     \         * ( QADV(IX,IZ+1,IDF) - QADV(IX,IZ,IDF) )
     \     - ( DENS0(IZ-1) + DENS0(IZ) ) / 2.0D0 
     \         * CDQNLZ(IX,IZ  ) 
     \         * ( QADV(IX,IZ,IDF) - QADV(IX,IZ-1,IDF) ) )
     \    / DENS0(IZ) / FDZ(IZ)

         DTQDNL(IX,IZ) = ( dqdinx + dqdinz ) / dx ** 2

   31 CONTINUE
C
      CALL BOUND( DTQDNL, NX, NZ, NXB, NZB, 1, 1, 1, 1 )
C
C================================================================
C
      RETURN
C
C
C
C===============================================================
      entry nldifi ( zgrid, nz )
C===============================================================
C
      if ( nz. gt. nzmax ) then 
         write(0,*) ' NLDIFA: nz too large '
         stop
      end if
      
      idturb = 1
      idtpot = 2
      idqvap = 3
      idqclw = 4
      idqrai = 5

*     for turbulent diffusivity 

      do iz = 0, nz-1
         sclap(iz,idturb) = 2000.0d0
         sclap(iz,idtpot) = 0.5d0
         sclap(iz,idqvap) = 0.01d0
         sclap(iz,idqclw) = 0.01d0
         sclap(iz,idqrai) = 0.01d0
      end do

c-------------------------------------------------------------------
      return



      E N D
